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ABSTRACT 


The  design  and  analysis  of  recoil  systems  for  direct  fire  weapons  has  been  conducted  at  Benet 
Labs  for  the  last  25  years.  The  model  developed  employs  the  Bernoulli  Equation  to  establish  the 
relationship  between  flow  rate  and  back  pressure  within  the  internal  paths  and  through  the  primary  orifice 
of  the  recoil  brake.  This  displacement  varying  orifice  provides  the  throttling  mechanism  needed  to 
generate  back  pressure  which  opposes  the  ballistic  driving  load  and  arrests  the  gun  in  recoil.  An  orifice  is 
designed  to  maintain  constant  upstream  pressure  over  the  complete  length  of  recoil,  thus  minimizing  the 
load  transferred  to  the  gun  support  and  vehicle.  The  equation  which  models  this  type  of  orifice  requires  an 
discharge  coefficient  (Cd)  which  ‘lumps’  all  of  the  flow  losses  due  to  contraction  and  directional  change  of 
the  fluid  stream.  For  the  model  developed  at  Benet  Labs,  this  coefficient  is  constant  regardless  of  fluid 
properties,  flow  regimes  and  geometries.  Test  data  from  firing  tests  and  research  experiments  indicate 
that  this  may  not  be  the  case. 

This  report  presents  the  details  of  using  an  offline  CFD  analysis  to  establish  the  flow  response 
characteristics  of  a  typical  recoil  brake  orifice  and  a  methodology  of  incorporating  these  results  via  a 
lookup  table  of  Cd  values  into  the  lumped  parameter  recoil  analysis  model  developed  at  Benet  Labs. 
Although  CFD  analysis  usually  requires  a  considerable  amount  of  time  and  cost  to  conduct,  we  have 
discovered  that  similarities  and  trajectories  in  the  solutions  allow  the  use  of  a  reduced  number  of  analysis 
conditions  between  flow  rates  and  geometries.  For  our  case,  this  type  of  CFD  analysis  should  become  a 
cost  competitive  method  for  use  in  hardware  design  and  system  analysis  when  the  conduction  of  a  large 
number  of  parametric  studies  is  required. 

Although  the  report  presented  herein  is  complete  in  terms  of  the  intent  of  the  research  proposal,  a 
considerable  amount  of  follow  on  studies  has  been  proposed  to  further  advance  the  use  of  these 
methodologies. 
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Introduction 


Recoil  Analysis  at  Benet  Labs:  Historical  Perspective  and  Modeling  Features 

In  the  early  1 980’s  Benet  Labs  was  given  the  mission  to  design  and  develop  tank  mounts  and 
recoil  systems  for  future  weapons  which  would  be  much  lighter  than  their  predecessors.  Several 
concepts,  computer  aided  routines  and  other  procedures  were  developed  to  fulfill  this  requirement.  The 
first  program  to  which  this  design  philosophy  was  applied  was  the  LW105mm  EX35  Weapon  (type 
classified  M35  in  1995).  Since  then,  there  have  been  several  developmental  weapons  to  which  Benet’s 
recoil  design  and  analysis  technique  has  been  applied.  Currently,  Benet  Labs  is  responsible  for  the 
design  of  the  primary  weapon  for  Future  Combat  System  (FCS).  It  is  a  1 20mm  smooth  bore  gun  tube,  a 
multi-lug  breech  and  a  box-type  mount  which  envelopes  two  recoil  brake  and  two  recuperator  assemblies. 
A  solid  body  representation  is  found  in  Figure  1.  Figure  2  presents  an  exploded  view  of  the  mount 
components.  The  two  green  components  are  the  recoil  brakes.  These  brakes  are  of  the  Schneider-type 
design  the  characteristics  of  which  may  be  found  in  [1]  and  [3]. 

The  function  of  a  recoil  system  is  to  transform  high  intensity  interior  ballistic  force  acting  in  the 
axial  direction  on  the  gun  tube  for  a  short  period  of  time  (5  to  20  ms),  into  much  lower  recoil  forces  acting 
on  the  support  structure  for  a  longer  period  of  time  (100  to  300  ms).  The  objective  is  the  attenuation  of 
short  duration  high  peak  loads  into  longer  duration  loads  whose  peak  values  are  much  lower.  In  addition, 
the  recoil  system  must  provide  for  the  storage  of  enough  energy  such  that  the  recoiling  parts  may  be 
returned  to  their  'home'  position  ready  for  the  next  firing.  Figure  3  contains  a  schematic  diagram  of  a 
weapon,  vehicle  and  the  recoil  system.  There  are  three  components  shown  in  the  upper  part  of  the  figure: 
1-  recoiling  parts;  2-  stationary  parts;  3-  recoil  system.  The  recoil  system  may  be  further  separated  into  a 
recoil  brake  and  recuperator.  Relative  force  distributions  in  time  which  drive  and  resist  are  shown  below 
the  components. 

The  function  of  the  recoil  brake  is  to  provide  the  bulk  of  the  retarding  forces.  This  is  accomplished 
by  metering  flow  through  a  travel  dependent  orifice.  The  metering  process  causes  upstream  pressure  to 
develop,  which  opposes  the  forces  on  and  momentum  of  the  recoiling  mass.  The  orifice  is  designed  such 
that  the  force  response  is  rather  constant  over  time  and  operating  stroke.  The  energy  absorbed  during 
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Figure  1 .  120mm  XM360  Primary  Weapon  Assembly 
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Figure  2.  Mount  Components 
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Figure  3.  Recoil  System  Schematic 


4 


this  process  is  converted  into  heat,  which  is  manifest  by  a  temperature  rise  in  the  brake  fluid  and 
structure.  A  schematic  of  a  generic  Schneider-type  recoil  brake  in  its  static  state  is  shown  on  Figure  4. 

The  main  components  of  the  assembly  are  called  out  on  the  figure.  Hydraulic  fluid  fills  the  internal 
volume.  As  the  brake  is  extended  due  to  ballistic  force  at  the  breech,  the  breech  end  cap,  cylinder,  control 
rod  and  buffer  valve  move  as  a  unit  in  the  recoil  direction.  Fluid  from  the  annular  chamber  between  the 
outer  diameter  (OD)  of  the  recoil  rod  and  the  inner  diameter  (ID)  of  the  cylinder  is  driven  along  the  OD  of 
the  recoil  rod  and  through  the  annular  opening  between  the  orifice  and  control  rod.  The  restricted  area 
between  the  orifice  and  control  rod  provides  the  impetus  for  the  development  of  the  upstream  pressure 
which  opposes  the  ballistic  driving  force  and  eventually  arrests  the  gun’s  recoiling  parts.  This  pressure  is 
a  function  of  this  resreictor  area  and  the  flow  rate  of  the  fluid.  As  shown,  the  control  rod  is  profiled  such 
that  the  pressure  distribution  in  this  chamber  remains  relatively  constant  as  the  velocity  of  the  recoiling 
parts  continually  change  as  the  gun  approaches  and  reaches  full  extension  in  recoil.  In  addition,  some  of 
the  fluid  flows  back  over  the  buffer  valve  into  the  counter  recoil  buffer  chamber  remaining  in  reserve  to 
decelerate  the  brake  as  it  returns  to  its  static  position.  Figure  5  contains  three  views  of  the  brake  during 
the  recoil  operation.  The  curved  arrows  indicate  the  flow  directions  from  the  main  annular  chamber, 
through  the  annular  area  between  the  control  rod  and  the  orifice  and  back  over  the  buffer  plug  and  into 
the  buffering  chamber  for  counter  recoil  control. 

Recoil  analysis  at  Benet  labs  is  accomplished  using  the  computer  code  RECOILA.  This  code 
written  in  FORTRAN  was  initially  developed  in  1982  and  has  continually  been  improved  through  several 
modifications.  It  employs  a  time  marching  integration  scheme  using  Adams  &  Bashfort  Four  Step  Method 
[4]  during  which  the  derivatives  from  three  preceding  time  steps  are  used  to  calculate  the  integral  solution 
at  the  current  time  step.  Using  the  applied  force  (ballistic  pressure)  and  the  resisting  forces  (brake 
pressure,  recuperator  pressure  and  friction)  the  equation  of  motion  is  solved  the  results  of  which  are 
reported  in  output  files  at  user  selected  time  steps.  The  heart  of  the  program  revolves  around  determining 
back  pressure  on  the  upstream  side  of  the  flow  within  the  brake  cylinder  to  determine  the  resistance  force 
applied  to  the  recoiling  parts.  This  quantity  is  computed  at  every  time  step  in  the  analysis.  Figure  6 
contains  a  schematic  of  flow  through  the  main  orifice  along  with  the  modeling  equations  used  to 
determine  the  relationship  between  flow  speed  and  back  pressure.  It  is  a  rather  simple  relationship 
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Figure  4.  Schneider  Recoil  Brake  Schematic  with  Components 
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derived  from  the  Bernoulli  equation  for  flow  through  a  pipe  with  restrictors.  Coupling  this  equation  with  the 
continuity  equation  for  incompressible  flow  yields  the  upstream  pressure  equation  (p^  which  is  a  function 
of  recoil  speed  (Vi),  area  ratio  (Ar  =  A-|/A2),  fluid  density  (p)  and  the  orifice  coefficient  (Cd).  The  current 
version  of  the  code  uses  a  constant  value  for  the  discharge  coefficient  (Cd)  based  solely  upon  textbook 
values  for  a  round-edge  orifice.  Published  values  for  this  configuration  are  about  0.95  [5].  It  is  the  feeling 
of  the  authors  that  this  parameter  is  a  strong  function  of  the  flow  geometry  as  well  as  fluid  properties  and 
kinematics  of  the  flow.  The  modeling  version  being  proposed  herein  will  have  the  capability  of  using  a 
variable  coefficient  for  flow  through  the  main  orifice  during  recoil.  We  expect  that  this  coefficient  will  be  a 
function  of  recoil  velocity  (Vi),  area  ratio  (Ar)  and  possibly  fluid  temperature. 

Current  Analysis  Techniques  Applied  to  XM360  Weapon 

Since  the  inception  of  the  recoil  design  and  analysis  mission  at  Benet  Labs,  thousands  of 
analysis  runs  have  been  conducted  and  compared  against  live  fire  test  data.  As  our  knowledge  regarding 
the  nuances  of  fluid  dynamics  within  damper-type  components  grew,  modifications  to  the  analysis 
algorithms  were  developed  as  well  as  additional  flow  features.  The  latest  version  (release  date:  October 
2006)  incorporates  fluid  compressibility,  where  the  density  is  treated  as  a  function  of  pressure  through  the 
bulk  modulus  of  the  fluid.  It  proved  to  be  a  minor  contributor  to  the  brake  pressure  response  (typically  a 
few  hundred  psi  for  a  nominal  pressure  of  4000  psi).  However,  the  use  of  a  constant  Cd  for  flow  through 
the  main  orifice  has  been  a  feature  of  this  code  since  the  first  release.  Discussion  of  a  dedicated 
laboratory  test  [6]  led  to  a  very  complicated  and  expensive  test  rig  which  may  or  may  not  have  generated 
the  data  needed  for  variable  Cd  confirmation.  This  avenue  was  abandoned  and  to  this  day  we  are  still 
using  a  constant  Cd  value  that  provides  acceptable  results  as  we  shall  subsequently  present  in  some 
detail. 


We  have  applied  the  current  analysis  methodology  to  test  data  generated  during  the  LW120mm 
XM360  ATD-4  Test  Program.  The  brake  pressure  results  (data  and  simulation)  for  Shot  #108  from  the  fall 
of  2006  are  shown  on  figures  7  and  8.  Figure  7  shows  pressure  data  (red  line)  and  simulation  response 
(blue  line)  for  both  brakes  plotted  against  time.  For  the  simulation,  the  measured  control  rod  profile  for 
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Figure  7.  Brake  Pressure  vs  Time  for  Round  #108  (Nominal  Cd  Value) 
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Figure  8.  Brake  Pressure  vs  Travel  for  Round  #108  (Nominal  Cd  Value) 
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each  brake  was  used  with  a  nominal  Cd  of  0.95  for  the  orifice.  For  both  brakes,  the  data  indicates  a  very 
similar  response.  The  peak  pressure  for  the  left  brake  is  slightly  greater  than  the  same  for  the  right  brake. 
Pressure  decay  after  peak  is  very  similar  in  distribution.  For  the  simulated  response,  slight  deviations  are 
shown.  Early  in  the  cycle  (0  to  10  ms)  the  simulation  under-predicts  the  data  however,  the  pressure  rise 
(10  to  15  ms)  is  quite  steep  achieving  a  sharp  peak  value  that  is  slightly  earlier  than  the  data  indicates. 
The  decrease  in  pressure  after  peak  is  very  quick.  This  is  most  likely  due  to  the  muzzle  brake  model  used 
in  the  simulation.  Muzzle  brake  parameters  are  simply  the  time  of  exit  of  the  projectile  and  the  Beta  value 
of  the  brake.  The  Beta  value  abruptly  decreases  the  gun’s  driving  load  (i.e.  ballistic  pressure)  subsequent 
to  exit.  So,  a  ‘step-like’  reduction  in  the  force  occurs  at  time  of  projectile  exit  resulting  in  an  abrupt  change 
in  recoil  dynamics.  After  peak,  the  response  settles  down  and  tracks  the  data  until  approximately  30  ms. 
Beyond  this  point,  the  simulation  deviates  significantly  from  the  data  meandering  above  and  below. 
Towards  the  end  of  the  recoil  cycle  the  simulation  response  indicates  pressure  values  that  are  several 
hundred  psi  greater  than  the  data  until  the  very  end  when  the  simulated  pressure  response  quickly  drops 
below  the  data  and  terminates  about  10  ms  earlier.  On  figure  8,  the  same  data  and  responses  are  plotted 
against  recoil  travel.  When  viewed  in  this  manner  the  deviations  between  the  data  and  simulation  are 
quite  clear.  From  0  to  2.5  inches,  the  simulated  response  is  less  than  the  data  by  300  to  500  psi.  The 
abrupt  rise  in  the  simulation  response  at  2.5  inches  is  very  clearly  shown.  At  its  peak  value,  the  simulation 
predicts  a  pressure  response  that  is  1000  to  1200  psi  greater  than  the  data  indicates  at  that  location. 

From  5  to  about  15  inches,  simulation  and  data  track  each  other.  From  15  to  23  inches,  the  simulated 
pressure  response  is  greater  than  the  data  by  200  to  500  psi. 

To  indicate  the  sensitivity  of  the  discharge  coefficient  we  conducted  a  series  of  simulations  of 
Shot  #108  using  various  values  for  Cd  and  compared  the  results  to  each  other  and  the  test  data.  Figures 
9  and  1 0  contain  graphical  representations  of  this  analysis.  We  have  used  four  values  for  Cd  (0.85,  0.90, 
0.95,  and  1 .00).  Pressure  plotted  against  time  is  shown  on  figure  9.  As  indicated,  for  low  values  of  Cd  the 
simulated  pressure  tracks  the  data  for  about  5  ms  then  rises  sharply  to  a  peak  which  surpasses  the  data 
by  about  1000  psi.  It  returns  and  tracks  the  data  from  40  to  60  ms  to  the  end  of  the  cycle.  For  large  values 
of  Cd  the  early  portion  of  the  cycle  is  matched  quite  well,  however,  towards  the  end  the  simulation  value 
surpasses  the  data  by  as  much  as  1200  psi.  The  distribution  in  time  from  60  to  1 10  ms  is  very  different  as 
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Figure  9.  Brake  Pressure  vs  Time  for  Round  #108  (Various  Cd  Values) 
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Figure  10.  Brake  Pressure  vs  Travel  for  Round  #108  (Various  Cd  Values) 
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well.  Pressure  plotted  against  travel  is  shown  on  figure  10.  The  differences  between  the  various  Cd  values 
and  the  data  are  easy  to  discern  when  viewed  in  this  manner.  It  is  quite  obvious  that  there  is  no  single 
value  for  Cd  that  best  matches  the  data.  Early  in  the  cycle  (0  to  2.5  inches)  a  Cd  value  of  0.85  is  evident. 
As  travel  progresses  (from  2.5  to  1 0  inches)  a  Cd  value  of  approximately  0.95  provides  the  best  fit. 
Towards  the  end  of  travel  (15  to  24  inches)  Cd  values  of  0.85  to  0.90  seem  to  best  match  the  data. 
Although  only  one  particular  shot  has  been  reported  herein,  this  trend  is  repeated  for  all  of  the  data  we 
have  simulated  to  date. 

Attempt  to  Extract  Discharge  Coefficients  from  Test  Data 

As  our  data  collection  methods  matured  especially  during  the  various  XM36  and  XM360 
programs,  we  have  been  able  to  get  a  much  better  and  more  accurate  reporting  of  the  recoil  system 
responses.  In  addition,  there  are  several  computer  programs  that  have  rendered  data  collection  and 
analysis  a  nearly  turnkey  operation.  There  is  one  caveat,  however,  data  analysis  and  review  should  be 
conducted  by  a  trained  professional  who  knows  the  expected  values  and  the  errors  associated  with  the 
measurement  transducers  and  reporting  hardware.  In  this  section,  we  attempt  to  reproduce  via 
assumptions  and  calculations  a  measure  of  the  instantaneous  Cd  value  by  using  test  data  from  the 
aforementioned  ATD-4  test. 

If  we  are  to  assume  that  fluid  flow  within  the  internal  chamber  of  the  brake  during  a  recoil  cycle 
follows  the  flow  laws  of  Bernoulli,  we  may  pick  any  two  points  in  the  flow  and  write  the  momentum  and 
continuity  equations  as  shown  in  figure  6  and  repeated  in  figure  11.  Note  that  Cd  is  buried  on  the  right 
side  of  the  flow  equation.  It  may  be  extracted  as  shown  in  the  final  equation  on  figure  1 1 .  Per  this 
equation,  Cd  is  directly  proportional  to  the  ratio  of  the  main  flow  area  (Ai)  to  the  area  at  the  throat  (A2), 
inversely  proportional  to  the  pressure  with  fluid  density  and  flow  speed  in  both  numerator  and 
denominator  within  the  radical.  Essentially,  knowing  the  pressure,  recoil  velocity  and  flow  areas  (orifice 
flow  area  is  a  function  of  travel)  we  may  in  theory  calculate  values  for  this  coefficient  using  the  time 
dependent  data  collected  during  the  recoil  portion  of  any  shot  in  a  firing  test.  We  have  routinely  collected 
pressure  and  travel  data  as  a  part  of  all  previous  XM36  and  XM360  firings,  however,  the  use  of  this  data 
to  calculate  Cd  is  problematic.  First,  since  a  parallel  flow  path  exists  to  the  buffer  pocket  we  do  not  know 
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Figure  11.  Extraction  of  Cd  from  Firing  Data 
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how  the  flow  is  apportioned  between  the  path  to  the  buffer  pocket  and  through  the  main  orifice.  Therefore, 
the  continuity  equation  on  the  figure  may  not  apply  as  shown.  A  potential  solution  to  this  problem  would 
be  to  ‘calculate’  flow  to  the  buffer  using  a  backflow  orifice  equation  (as  a  flow  restriction  along  this  path) 
along  with  the  pressure  data.  Once  this  flow  rate  to  the  buffer  pocket  is  estimated  the  remaining  fluid  must 
pass  through  the  main  orifice.  Second,  although  the  travel  is  globally  representative  of  the  gun’s  position 
it  is  not  smooth  enough  to  yield  an  accurate  recoil  velocity  using  a  point  by  point  differentiation  algorithm. 
A  potential  solution  to  this  problem  would  be  to  fit  the  travel  data  to  a  function  using  a  least  squares 
method  and  then  differentiate  the  function  to  determine  velocity.  Both  MatLab®  and  Sigma  Plot®  have 
several  standard  functions  to  use  in  a  least  squares  fit  algorithm.  After  several  tries  of  both  algebraic  and 
transcendental  functions  the  existence  of  3  sigmoidal-type  functions  within  the  Sigma  Plot  regression 
library  seems  to  best  fit  the  data  and  expected  velocity  profiles. 

The  first  is  the  3  parameter  Logistic  function,  which  is  depicted  on  Figure  12  for  Shot  #104.  This 
figure  contains  2  plots  and  the  relevant  fitting  equations.  The  upper  plot  contains  travel  data  from  the  test 
with  the  logistic  fitting  function  superimposed.  On  the  lower  plot,  the  velocity  as  derived  from  an  analytic 
differentiation  of  the  function  is  shown.  As  indicated  in  the  equation,  the  travel  (x)  is  a  function  of  time 
which  appears  only  in  the  denominator  of  the  function  along  with  the  parameters  a,  b  and  t0.  The 
residuals  (the  difference  between  the  data  and  the  fit  function)  are  on  the  order  of  ±0.2  inches.  The 
meandering  nature  of  the  data  is  observed  through  the  slight  oscillation  about  the  smoother  fit  function. 

As  mentioned  earlier,  this  is  most  likely  due  to  the  transducer  used  for  travel  determination  and  is  not 
unlike  other  data  that  has  been  collected  before.  In  the  velocity  function  ( v )  shown  on  in  the  lower  plot, 
the  independent  variable  (#)  appears  in  both  the  numerator  and  denominator  along  with  combinations  of 
the  other  parameters.  We  expect  a  velocity  response  of  this  type;  however,  the  actual  shape  of  the  tail 
end  (i.e.  the  decay  portion)  is  unknown.  The  second  is  the  3  parameter  Chapman  function,  which  is 
depicted  on  Figure  13  for  Shot  #104.  As  indicated  in  the  equation,  the  travel  (x)  is  a  strong  function  of  the 
exponential  (e).  Since  the  exponential  appears  in  both  the  numerator  and  denominator  of  the  velocity 
function,  its  value  at  zero  time  is  indeterminate;  however,  by  using  techniques  of  calculus,  it  is  found  to  be 
zero.  The  third  is  the  3  parameter  Hill  function,  which  is  depicted  on  Figure  14  for  Shot  #104.  As  indicated 
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Figure  12.  Data  and  Logistic  Fit  of  Recoil  Travel  for  ATD-4  Shot  #104 
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Figure  13.  Data  and  Chapman  Fit  of  Recoil  Travel  for  ATD-4  Shot  #104 
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Figure  14.  Data  and  Hill  Fit  of  Recoil  Travel  for  ATD-4  Shot  #104 
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in  the  equation  the  travel  (x)  and  velocity  (v)  are  rational  functions  of  (f).  The  form  of  the  velocity  equation 
admits  a  value  of  zero  at  zero  time.  In  summary,  all  three  fitting  methods  produced  travel  and  velocity 
data  that  are  nearly  identical.  We  shall  proceed  using  all  three  methods  to  analyze  the  firing  data  in  the  Cd 
equation  from  figure  1 1 . 

A  MatLab®  script  (see  Appendix  A)  was  developed  and  used  for  processing  the  recoil  pressure, 
brake  geometry,  travel  and  velocity  fitting  functions  to  generate  the  dynamic  Cd  values.  Each  data  set 
contained  about  1100  points  for  each  variable.  We  shall  discuss  in  some  the  detail  the  results  for  this 
processing  of  Round  #104  and  highlight  the  remaining  10  rounds  in  this  series.  As  mentioned  earlier,  in 
order  to  use  firing  data  for  this  process,  it  must  be  rather  smooth  and  continuous.  As  shown,  the  raw 
travel  data  did  not  possess  this  characteristic;  therefore,  a  smoothing  technique  had  to  be  implemented  to 
determine  recoil  velocity.  Additionally,  the  brake  pressure  data  which  comprises  a  term  in  the  Cd 
calculation  was  rather  oscillatory  as  well,  especially  early  in  the  cycle,  therefore  a  smoothing  process  had 
to  be  applied.  Since  we  only  needed  the  data  and  not  its  derivative,  a  polynomial  function  would  be  good 
enough  for  its  representation.  On  Figure  15,  the  raw  data  and  its  15  term  polynomial  fitting  function  is 
shown  for  Round  #104.  The  red  line  is  the  data  and  as  shown  it  contains  low  amplitude  oscillations  up  to 
about  35  ms.  The  fitting  function  (blue  line)  yields  a  very  good  representation  of  the  data,  smoothing  the 
initial  transients  very  well. 

On  Figure  1 6,  the  result  for  the  Cd  calculation  is  shown  for  the  Logistic  fit  of  the  travel  data.  This 
coefficient  is  plotted  against  recoil  travel  using  the  pressure  data  for  the  left  and  right  brakes 
independently.  The  upper  plot  contains  the  results  for  left  brake  while  the  result  for  the  right  brake  is 
shown  on  the  lower  plot.  The  result  for  either  brake  is  nearly  the  same.  Up  to  about  2.5  inches  of  travel 
the  Cd  value  is  low  (-0.60  to  0.70).  During  this  portion  of  travel,  the  speed  is  slowly  rising  and  the  throat 
area  is  constant,  however,  Cd  rises  quickly,  reaching  a  maximum  value  of  slightly  above  1.1  at  5  inches  of 
travel.  Obviously  a  value  of  1.1  is  not  physically  possible  since  Cd  values  should  be  between  0  and  1. 
However,  during  this  portion  of  recoil,  the  speed  and  pressure  are  increasing  rapidly,  whereas,  the  throat 
area  is  decreasing,  so  the  calculated  value  (since  a  straightforward  equation  is  being  used)  could  rise 
above  1.0  although  this  is  ambiguous.  Maximum  recoil  velocity  and  pressure  are  shown  to  occur  at  about 
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Figure  15.  Data  &  Polynomial  Fit  of  Brake  Pressure  for  ATD-4  Shot  #104 
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CALCULATED  Cd  from  RIGHT  BRAKE  vs  TRAVEL 


Figure  16.  Cd  Calculation  Using  Logistic  Data  Fit  for  ATD-4  Shot  #104 
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5  inches.  The  Cd  value  reverses  after  the  5  inch  location  reducing  to  values  slightly  less  than  1 .0  up  to 
about  the  10  inch  travel  location.  From  10  to  about  20  inches  of  travel,  the  Cd  value  monotonically  drops 
to  less  than  0.8.  Beyond  the  20  inch  location,  the  value  rises  rapidly  to  values  in  excess  of  1 .2.  During  this 
phase  of  the  cycle,  pressure,  velocity  and  throat  area  are  decreasing  quickly,  therefore,  the  transient 
nature  of  the  variables  involved  in  the  calculation  may  render  its  result  meaningless.  Bottom  line,  once 
flow  is  established  and  the  variables  are  changing  slowly  (from  5  to  20  inches  of  travel)  the  calculation  of 
Cd  is  probably  quite  good.  Nearly  the  same  results  and  comments  may  be  made  for  the  Chapman  and  Hill 
fit  which  are  shown  on  Figures  17  and  18.  This  similarity  in  response  carries  throughout  the  remaining  10 
rounds  of  this  test.  This  analysis  lends  credence  to  the  premise  that  Cd  is  variable  and  probably  a  function 
of  flow  speed  and  the  ratio  of  upstream  chamber  area  to  the  throat  area. 

Use  of  Alternate  Methodologies  for  Calculating  Discharge  Coefficient  Values 

The  idea  that  Cd  values  are  dependent  upon  various  fluid  and  geometry  parameters  for  flow 
around  the  circumference  of  plates  mounted  within  a  cylinder  was  shown  experimentally  by  Bell  and 
Bergelin  [7].  The  basic  test  fixture  was  a  long  cylinder  with  an  inner  diameter  of  5.2541  inches  within 
which  individual  plates  of  sizes  between  5.0251  to  5.2325  inches  in  diameter  of  varying  edge 
configurations  and  thickness  were  mounted  both  concentrically  and  tangentially.  Three  configurations  of 
the  plates’  outer  circumference  were  tested.  They  were  square  edge  with  measurable  flow  length,  sharp 
edge  and  rounded  entrance  edge.  Both  oil  and  water  were  used  and  as  it  was  pumped  through  the  test 
section  (after  thermal  and  flow  equilibrium  was  reached),  pressure  and  flow  rate  data  was  collected.  The 
orifice  coefficients  for  each  trial  was  calculated  based  upon  the  test  data  and  a  standard  equation  from  a 
Bernoulli  flow  analysis  similar  to  that  explained  in  the  last  section.  Additionally,  the  Reynolds  number  of  all 
experimental  trials  was  calculated  using  orifice  and  cylinder  diameters  as  the  length  parameter  along  with 
the  fluid’s  density,  viscosity  and  flow  rate. 

The  results  are  presented  in  the  form  of  log-log  plots  of  Cd  as  a  function  of  Reynolds  Number. 
Basically  two  distinct  responses  developed;  one  for  the  sharp  edge  orifice  and  another  for  the  square  and 
round  edge  orifice.  These  are  presented  in  two  forms  on  figures  5  through  figure  8  of  the  referenced 
report.  In  addition,  several  empirical  equations  relating  the  test  parameters  (i.e.  orifice  geometry  and 
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LW120  XM36  ATD-4  RECOIL  SYSTEM  TEST  RESPONSE 
TEST  ROUND  #104  (3  PARAMETER  CHAPMAN  FIT) 


CALCULATED  Cd  from  LEFT  BRAKE  vs  TRAVEL 


CALCULATED  Cd  from  RIGHT  BRAKE  vs  TRAVEL 


Figure  17.  Cd  Calculation  Using  Chapman  Data  Fit  for  ATD-4  Shot  #104 
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DISCHARGE  COEFF.  (-)  DISCHARGE  COEFF. 


LW120  XM36  ATD-4  RECOIL  SYSTEM  TEST  RESPONSE 
TEST  ROUND  #104  (3  PARAMETER  HILL  FIT) 


CALCULATED  Cd  from  LEFT  BRAKE  vs  TRAVEL 


CALCULATED  Cd  from  RIGHT  BRAKE  vs  TRAVEL 


Figure  18.  Cd  Calculation  Using  Hill  Data  Fit  for  ATD-4  Shot  #104 
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Reynolds  number)  to  Cd  were  developed,  the  graphs  of  which  were  superimposed  on  the 
abovementioned  plots.  To  summarize,  Cd  tends  to  increase  monotonically  with  increasing  Reynolds 
Number  (within  a  2  to  20,000  range)  and  asymptotically  approaching  a  limiting  value  dependent  upon  the 
edge  geometry  of  the  disk. 

Although  this  work  is  of  great  value  for  the  geometries  that  were  tested  we  could  not  apply  these 
results  to  the  problem  at  hand.  Our  flow  restrictor  was  much  smaller  in  size  than  that  of  the  test  and 
contains  a  moving  part  (Control  Rod)  which  has  a  variable  outer  diameter  and  moves  through  the  orifice 
at  a  varying  speed.  Additionally,  the  orifice  is  fed  from  5  ports  in  the  piston  head  each  of  which  has  a  90° 
shoulder  of  0.31  inch  at  the  leading  edge  of  the  entrance  and  a  diverter  angle  of  approximately  50°.  In  the 
past,  discussions  regarding  laboratory  testing  of  the  orifice  configuration  was  discussed  [6],  but  as 
mentioned  earlier,  the  parameters  needed  to  conduct  experiments  in  the  flow  regime  of  a  live  fire  test 
could  not  be  achieved.  Additionally,  a  first  cut  at  an  experimental  setup  seemed  to  be  quite  pricey  and 
was  eventually  abandoned. 

In  the  late  1990’s,  Benet  Labs  began  using  dedicated  computational  fluid  dynamics  (CFD)  codes 
for  analysis  of  gas  flow  through  muzzle  brakes  to  predict  overpressure  fields  in  close  proximity  to  the 
muzzle  end  of  the  gun.  These  fields  were  the  result  of  expulsion  and  expansion  of  the  propellant  gases 
after  the  projectile  has  exited  the  muzzle.  This  type  of  flow  is  highly  complicated  in  that  gas  compression, 
shock  wave  development  and  composition  changes  occur  in  a  highly  transient  environment.  These 
analyses  were  applied  to  several  gun  systems,  including  the  XM360.  Informal  results  indicate  that  the 
correlation  between  analysis  and  test  data  is  good  (although  costly  in  terms  of  execution  times).  In 
addition,  CFD  was  applied  to  an  unsteady  simplified  2-D  axisymmetric  model  of  a  recoil  cylinder  buffer 
plug,  control  rod,  and  orifice  inlet  [8].  This  model  utilized  dynamic  meshing  to  simulate  motion  of  the 
recoil  cylinder.  This  unsteady  model  utilized  oil  as  a  medium,  but  also  simulated  the  multi-phase  flow 
environment  to  account  for  cavitation  and  collapse  of  the  fluid.  Based  on  these  results,  it  makes  sense 
that  the  application  of  CFD  modeling  for  an  incompressible  Newtonian  fluid  flow  in  the  absence  of  heat 
transfer,  phase  changes,  etc.  should  be  relatively  low  cost.  Additionally,  if  a  representative  axisymmetric 
version  of  the  flow  path  could  be  conceived  this  would  render  the  problem  to  be  two-dimensional,  further 
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reducing  run  times.  The  combination  of  these  three  distinct  bodies  of  work  cited  above  was  the 
motivation  to  conduct  detailed  CFD  analysis  of  fluid  flow  within  the  brake  during  the  recoil  cycle.  We 
would  attempt  to  conduct  ‘numerical  experiments’  in  a  manner  similar  to  the  experimental  work  of  Bell  and 
Bergelin.  The  result  would  be  a  table  of  Cd  values  as  a  function  of  one  or  several  parameters  of  the  flow. 

Incompressible  Flow  Modeling  Using  CFD  Methodology 

© 

Description  of  CFD  in  General  and  the  Code  CFdesiqn 

Finite  element  analysis  methods  (fea)  have  been  used  to  analyze  elastic  structures  for  several 
decades.  In  general,  fea  models  partition  the  analysis  domain  into  many  ‘nodes’  at  which  displacements 
are  calculated  and  reported.  The  areas  or  volumes  between  the  nodes  are  filled  in  with  elements 
throughout  which  these  nodal  displacements  are  approximated  by  basis  functions  using  the  values  at 
adjacent  nodes  for  their  ‘data’.  The  displacements  vector  and  element  formulation  are  assembled  into  a 
matrix  expression  which  is  equated  to  a  vector  of  know  forces  much  like  a  standard  system  of  algebraic 
equations.  The  difference  is  that  the  solution  vector  (nodes)  is  quite  long  whereas  the  matrix  representing 
the  element  formulation  is  the  square  of  the  length  of  the  nodal  vector.  Solution  to  a  problem  formulated 
in  this  manner  is  accomplished  through  the  use  of  clever  numerical  techniques.  Once  the  displacement 
vector  has  been  solved,  the  basis  functions  and  the  known  values  for  the  mechanical  properties  of  the 
material  are  used  to  determine  strains  and  stresses  at  any  location  in  the  structure.  This  is  a  field  variable 
problem  where  the  stresses  are  the  field  variables.  One  may  liken  this  to  a  series  of  rings  attached  to 
each  other  by  springs.  If  one  or  many  of  the  rings  are  displaced,  the  springs  will  either  extend  or 
compress.  The  extension  or  compression  will  determine  the  internal  force  within  each  spring.  The  elastic 
energy  in  each  spring  is  directly  proportional  to  these  forces.  Stress  and  strain  is  directly  proportional  to 
these  forces.  For  the  types  of  structures  and  loads  imparted  to  components  designed  by  engineers  and 
technicians  at  Benet  Labs,  the  use  of  fea  analysis  is  crucial  to  the  success  of  weapons’  systems.  We 
have  been  using  fea  methods  for  design  for  nearly  4  decades.  The  initial  code  used  was  NASTRAN 
developed  by  NASA  for  the  aerospace  industry.  During  the  early  1980’s,  we  acquired  the  more 
comprehensive  and  adaptable  code  called  ABAQUS  which  contains  more  features  than  NASTRAN.  Also, 
there  are  several  other  specialty  fea  codes  that  have  been  developed  in-house  or  acquired  commercially 
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for  desktop  use. 


Finite  element  analysis  methods  are  not  limited  to  structural  problems  only.  A  growing  field  of 
interest  is  in  computational  fluid  dynamics  (CFD)  in  which  the  methodology  of  finite  elements  has  been 
incorporated  into  finite-volume  problems.  CFD  uses  setup  procedures  that  are  similar  to  the  structural 
models.  The  fluid  contained  in  a  spatial  domain  is  discretized  into  small  cells  (elements)  the  intersections 
of  which  are  the  nodes.  A  suitable  algorithm  to  solve  the  Navier-Stokes  equations  of  motion  is  then 
applied.  These  equations  describe  the  changes  in  momentum  within  an  infinitesimal  volume  of  fluid  as 
simply  the  product  of  changes  in  pressure  and  dissipative  viscous  forces  acting  inside  the  fluid.  These 
equations  are  differential  equations  which  describe  the  relationships  among  the  rates  of  change  of  the 
variables. 

CFdesign  is  one  such  CFD  code  that  can  be  run  locally  on  a  high  powered  PC.  Its  analysis 
method  is  based  upon  the  finite-element  method,  unlike  other  fluid  flow  analysis  packages  that  utilize  a 
finite-volume  approach.  This  method  lends  itself  to  a  more  flexible  representation  of  the  geometry  and  the 
physical  properties  of  the  fluid  and  components  in  the  analysis.  Additionally,  it  derives  the  geometric  data 
from  several  commercial  CAD  packages  (at  Benet  Labs  Pro/Engineer®  is  used).  Regarding  flow 
modeling,  it  has  the  ability  to  simulate  transient  and  steady  state,  compressible  or  incompressible,  laminar 
or  turbulent,  internal  or  external,  and  relative  motion  of  walls  [9].  For  the  problem  at  hand  (i.e.  pipe  flow 
within  a  recoil  brake)  we  chose  a  steady  state,  incompressible  analysis  with  the  possibility  of  both  laminar 
and  turbulent  zones. 

Application  to  the  Current  Problem  via  Simplified  Axisymmetric  Geometry 

The  current  Benet  Labs  recoil  analysis  code  uses  a  lumped  flow  resistor  to  model  the  effects  of 
the  shoulder,  angular  paths  and  restricted  area  of  the  orifice  with  no  requirement  to  calculate  the  detailed 
pressure  or  flow  distributions  throughout  these  volumes.  On  Figure  1 9,  the  assembly  view  of  a  recoil 
brake  is  shown.  The  cylinder  and  control  rod  move  to  the  right  driven  by  the  pressure  in  the  breech  which 
accelerates  the  gun  in  recoil.  The  orifice  which  is  attached  to  the  throttling  sleeve  and  piston  head  remain 
stationary.  Fluid  is  forced  through  the  annular  area  between  the  control  rod  profile  and  orifice.  The  path 
through  the  6  feed  ports  in  the  piston  head  and  the  restricted  area  between  the  orifice  and  control  rod  are 
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BRAKE  CYLINDER  CONTROL  ROD 


Figure  19.  Assembly  View  of  XM360  Recoil  Brake 
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the  areas  of  concern.  As  such,  the  actual  flow  is  forced  around  shoulders  and  diverted  through  angular 
paths  as  well  as  a  possible  swirling  motion  as  the  fluid  is  driven  circumferentially  towards  one  of  the  ports. 

Although  we  have  always  treated  this  as  axisymmetric  pipe  flow,  the  exclusion  of  the  3D 
component  through  the  six  ports  in  the  piston  head  should  not  compromise  the  required  accuracy  of  the 
model.  Therefore,  flow  representation  for  our  CFD  analysis  will  be  axisymmetric. 

A  simplified  model  of  the  process  detailed  above  is  shown  on  Figure  20.  It  is  comprised  of  two 
solid  boundaries  (heavy  dark  blue  line)  which  envelopes  the  fluid  (light  blue  fill).  The  model  is 
axisymmetric  and  the  given  radial  dimensions  best  represent  the  actual  specifications  for  the  current 
recoil  brake.  The  variable  dimension  is  the  rod  diameter,  shown  on  the  right  side  of  the  lower  boundary.  A 
constant  rod  diameter  was  utilized  rather  than  a  taper  to  simulate  the  control  rod.  This  constant  diameter 
varies  from  0.375  to  0.5625  inches.  We  choose  increments  of  0.01 0  inches  (ending  at  0.555  in  due  to  run 
time  constraints  as  will  be  detailed  later)  yielding  1 9  distinct  geometric  models.  The  length  of  the  model  is 
only  4.75  inches  which  encompasses  a  small  portion  of  the  total  flow  within  the  brake.  (When  fully 
compressed,  the  actual  length  of  the  fluid  column  in  the  high  pressure  volume  is  24  inches.)  Unlike  the 
physical  component,  the  solid  boundaries  are  immobile  whereas  the  fluid  is  fed  from  left  to  right  within  the 
enclosure.  The  fluid  entering  the  left  boundary  is  given  a  flow  rate  or  velocity  which  simulates  a  given 
recoil  speed  whereas  the  fluid  exiting  at  the  right  boundary  is  in  a  zero  pressure  state  which  is  the 
assumed  state  in  the  recoil  simulation  model.  The  CFD  model  calculates  back  pressure  near  the 
entrance.  Using  back  pressure  as  the  first  independent  variable  and  the  input  flow  speed  as  the  second 
independent  variable  the  reformulated  Bernoulli  equation  (Figure  1 1)  is  solved  for  Cd .  We  choose  speeds 
from  5  to  700  in/sec  in  varying  increments  as  our  inlet  flow  condition.  The  full  matrix  of  geometries  and 
flow  speeds  yielded  627  runs  of  the  CFD  code. 

The  element  density  is  a  function  of  the  location  within  the  fluid  and  the  clearance  between  the 
boundaries.  Prior  to  developing  this  model,  we  conducted  convergence  studies  on  a  similar  flow  regime 
and  concluded  that  a  minimum  of  8  triangular  elements  within  the  core  of  flow  and  3  prismatic  boundary 
layer  elements  at  each  surface  was  sufficient  to  produce  convergent  results  as  a  function  of  element 
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CENTERLINE  OF  CROSS  SECTION 
(ALL  VERTICAL  DIMENSIONS  ARE  RADIAL) 
(SOME  DIMENSIONS  OMITTED  FOR  CLARITY) 


Figure  20.  CFD  Model  of  Axisymmetric  Flow  Area 
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Figure  21 .  Zones  of  the  CFD  Model  with  Typical  Meshing  Density 


density.  Given  that  clearance  values  would  be  as  small  as  0.010  inches  for  a  tight  fitting  control  rod,  we 
had  the  model  divided  into  5  zones  within  which  the  element  density  could  be  independently  specified. 
These  zones  are  depicted  on  Figure  21  as  Z-1  through  Z-5.  The  zones  Z-1  and  Z-5  at  the  extreme  end  of 
the  fluid  volume  are,  for  the  most  part,  independent  of  throat  clearance;  therefore,  element  density  within 
these  two  zones  was  is  always  set  to  about  20  triangular  and  3  each  rectangular  boundary  layer  elements 
at  the  solid  surfaces.  Nodal  spacing  was  between  0.035  to  0.040  inches.  For  zones  Z-2  and  Z-4,  which 
must  merge  with  the  elements  within  the  throat  area,  there  is  some  dependence  upon  clearance.  Nodal 
spacing  was  0.0175  to  0.025  inches  dependent  upon  the  throat  clearance.  Spacing  in  this  range  yielded 
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element  densities  of  32  to  40  triangular  elements  along  with  the  usual  boundary  layer  elements.  Within 
the  throat  zone  (Z-5),  nodal  spacing  varied  from  0.0025  to  0.015  inches,  dependent  upon  clearance.  This 
range  always  kept  the  element  number  to  at  least  8  triangular  elements,  the  condition  needed  for 
convergence.  For  the  looser  fitting  control  rod  diameters,  the  number  of  fluid  elements  was  approximately 
20,000  whereas  for  the  tightest  fit  the  number  of  elements  increased  to  234,000.  For  the  mesh  shown  in 
the  figure,  there  are  11,349  nodes  and  21,970  elements.  The  nodal  density  change  from  Zones  1  to  3  is 
quite  evident.  Several  hundred  iterations  were  needed  for  convergence  irrespective  of  number  of 
elements.  The  actual  run  times  for  a  single  analysis  ranged  from  10  minutes  for  coarser  element  models 
to  over  8  hours  for  the  denser  element  models.  Several  weeks  were  needed  to  complete  the  full  range  of 
simulations. 

Representative  Results  from  the  CFD  Steady  State  Flow  Model 

As  mentioned,  we  conducted  approximately  627  individual  runs  to  generate  the  pressure  and  flow 
speed  data  needed  to  calculate  Cd.  We  shall  report  in  detail  the  results  from  one  geometric  model.  The 
chosen  configuration  is  the  one  in  which  the  rod  diameter  is  1.010  inches  (0.505  inches  radial).  This 
geometric  configuration  needed  a  total  of  43,450  elements  for  convergence.  At  this  element  density,  the 
compilation  time  per  iteration  was  2.6  seconds.  For  the  33  runs  the  total  number  of  iterations  was  1 8,950 
which  convert  to  about  14  hours  of  run  time.  Since  batch  processing  and  overnight  runs  are  possible,  all 
33  including  post  processing  were  completed  in  a  few  days.  All  the  numerical  data  required  to  determine 
Cd  was  ported  to  a  Sigma  Plot  worksheet  and  processed  via  a  user  developed  transform.  (A  transform  is 
merely  a  set  of  instructions  for  manipulating  and  processing  data  that  reside  in  Sigma  Plot  spreadsheet 
files.) 


On  Figure  22,  the  calculated  values  for  Cd  are  shown  in  two  plots.  The  upper  reports  Cd  as  a 
function  of  flow  velocity  which  is  the  inlet  condition  used  in  the  CFD  model,  whereas  the  lower  plot  shows 
Cd  as  a  function  of  volumetric  flow  rate.  (Velocity  and  flow-rate  at  any  point  in  the  fluid  are  directly 
proportional  via  the  continuity  equation  for  incompressible  flow.)  These  plots  are  typical  in  shape  for  all 
rod  profiles;  however,  the  Cd  values  are  shifted  vertically  as  a  function  of  rod  diameter.  From  these  plots  it 
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CFD  ANALYSIS  OF  RECOIL  CYLINDER  /  ORIFICE  FLOW 
FLUID:  INCOMPRESSIBLE  HYDRAULIC  BRAKE  @  75°F 


FLOW  RATE  (in3 /sec) 

Figure  22.  Discharge  Coefficient  vs  Velocity  &  Flow  Rate 
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is  apparent  that  Cd  is  a  strong  function  of  flow  speed.  At  low  flow  speeds  (0  to  1 00  in/sec)  the  rate  of 
change  is  quite  large  from  0.65  to  0.82.  From  this  point  it  rises  at  a  much  slower  rate,  finally  achieving  its 
highest  value  of  0.925  at  a  flow  speed  of  700  in/sec.  It  is,  as  expected,  a  well  behaved  response,  void  of 
singularities  and  inflection  points.  The  data  needed  for  the  recoil  analysis  code  will  be  derived  from  the 
lower  plot.  We  expect  to  set  up  tables  of  Cd  versus  flow  rate,  since  this  is  the  value  that  is  iterated  in  the 
recoil  model  when  attempting  to  calculated  correct  distribution  of  flow  rate  and  pressure. 

Figures  23  and  24  contain  more  detailed  information  regarding  the  flow.  On  Figure  23,  we  show 
the  axial  velocity  profile  from  the  entrance  to  the  exit  of  the  orifice.  The  range  of  values  as  shown  on  the 
color  bar  is  from  0  to  1 000  in/sec  with  the  value  of  700-1000  within  the  fluid’s  core  at  the  throat.  The  flow 
speed  at  the  inlet  is  only  25  in/sec,  so  a  speedup  of  nearly  40  to  1  is  indicated  based  upon  the  maximum 
flow  speed  in  the  throat.  By  using  the  continuity  relationship,  the  average  speed  increase  should  be  about 
30  to  1  for  this  orifice  which  is  the  case  if  we  were  to  calculate  an  average  velocity  at  the  throat.  On  figure 
24,  we  show  the  pressure  drop  through  the  orifice.  Left  of  the  shoulder  the  pressure  is  rather  constant 
between  30  and  35  psi.  Pressure  drop  does  not  occur  until  the  flow  becomes  quite  restricted  just  to  the 
left  of  the  throat,  finally  exiting  on  the  right  at  about  0  psi.  The  data  on  these  two  plots  represent  a 
classical  result  for  orifice  flow,  the  rapid  speed  up  and  pressure  drop  as  the  flow  path  becomes  more 
restricted. 

In  the  following  section,  we  shall  discuss  the  details  involved  in  compiling  these  results  into  a  form 
which  in  amendable  to  incorporation  into  the  Benet  Labs  recoil  analysis  code.  Additionally,  we  shall  delve 
into  the  logic  used  in  the  code  to  implement  a  variable  Cd  for  the  main  orifice. 
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Figure  23.  Velocity  Profile  Through  Orifice  505  @  25  in/sec  Boundary 
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(5)  Sialic  Pressure  -psi 


Figure  24.  Static  Pressure  Profile  Through  Orifice  505  @  25  in/sec  Boundary 


38 


Application  of  Results  to  the  Benet  Labs  Recoil  Analysis  Model 

Extraction  of  Discharge  Coefficient  Values  from  CFD  Results 

We  have  conducted  627  independent  runs  using  CFdesign®  code  for  19  geometry  configurations 
and  33  flow  speed  conditions.  Total  CPU  time  for  the  conduction  of  these  analyses  was  485  hours  with  a 
total  number  of  iterations  of  approximately  320,000.  In  the  last  section,  we  discussed  in  some  detail  the 
‘raw’  data  results  from  the  model  as  well  as  the  method  of  determining  the  parameter  of  concern  which  is 
the  discharge  coefficient  of  the  restrictive  flow  path.  Reiterating,  the  Benet  Labs  recoil  analysis  code  uses 
the  Bernoulli  equation  for  all  its  pressure  versus  flow  rate  calculations,  therefore,  any  restrictive  areas 
must  be  rated  for  their  efficiency  (orifice  coefficients)  with  respect  to  pressure  enhancement  due  to 
contraction  of  the  flow  stream.  Handbooks  have  a  plethora  of  data  regarding  these  coefficients  with  such 
parameter  features  as  the  lead  in  and  exit  angle  of  the  constriction.  The  orifice  type  used  in  brakes 
designed  at  Benet  Labs  utilizes  a  spherical  boundary  on  the  orifice  surface  and  a  gentle  taper  on  the 
control  rod  profile. 

From  the  results  cited  in  the  previous  section,  we  decided  that  the  best  way  to  attack  the  problem 
was  to  develop  a  two-parameter  dataset  to  represent  the  Cd  values.  The  parameters  chosen  are  the  flow 
rate  in  cubic  inches  per  second  (in3/sec)  and  a  parameter  called  the  area  ratio.  The  range  of  flow  rate 
values  used  in  the  CFD  analysis  was  31  to  4340  (in3/sec)  which  bounds  the  expected  recoil  velocities. 
The  area  ratio  is  the  ratio  of  the  flow  area  in  the  main  chamber  to  the  area  across  the  throat  of  the 
restrictor.  The  range  of  this  value  in  the  actual  brake  is  1 1  to  702.  When  the  area  ratio  is  at  a  low  value 
the  throat  area  is  open  and  vice  versa  when  the  ratio  is  large.  Results  for  7  of  the  19  separate  area  ratio 
models  are  shown  on  Figures  25  and  26  using  color  coded  lines  to  indicate  the  area  ratio.  The  ‘data’ 
presented  is  the  calculated  Cd  values  as  a  function  of  flow  rate  with  the  various  area  ratio  values  as 
parameters  of  each  trajectory.  Four  flow  rate  zones  were  established  based  on  the  corresponding  recoil 
velocities.  In  Zone  1  the  flow  rate  range  is  0  to  400  in3/sec.  For  Zone  2  the  flow  rate  is  400  to  1000 
in3/sec.  Results  for  these  two  flow  rates  are  shown  on  Figure  25.  Zone  3  flow  rate  is  1000  to  2600  in3/sec. 
whereas  Zone  4  flow  rate  is  2600  to  4400  in3/sec.  Results  for  these  rates  are  shown  on  Figure  26. 

The  upper  plot  on  Figure  25  contains  results  for  Zone  1  flow  rates.  The  range  of  flow  rate  values 
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Figure  25.  Cd  Results  for  Flow  Rate  Zones  1  &  2 
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Figure  26.  Cd  Results  for  Flow  Rate  Zones  3  &  4 
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is  0  to  400  in3/sec  which  corresponds  to  a  recoil  velocity  range  of  5  to  70  in/sec.  These  flow  rates  would 
be  achieved  very  early  in  the  recoil  cycle  when  the  initial  acceleration  causes  a  ramp  up  in  speed  and 
again  at  the  termination  of  the  cycle  when  deceleration  reduces  recoil  speed.  The  smaller  area  ratios 
would  be  the  case  in  the  early  portion  of  the  cycle  and  the  larger  for  the  final  portion  of  recoil  cycle.  The 
Cd  range  for  these  conditions  is  0.40  to  0.86.  As  indicated,  the  data  trajectories  are  smooth  and 
monotonic  with  respect  to  both  flow  rate  and  area  ratio.  A  substantial  ‘ramp  up’  is  shown  for  very  low  flow 
rates  (0  to  150  in3/sec).  The  Cd  trajectories  decrease  with  increasing  area  ratio  with  a  several  point 
decrease  for  the  largest  area  ratio  but  remain  rather  typical  in  shape.  The  lower  plot  on  Figure  25 
contains  results  for  Zone  2  flow  rates.  The  range  of  flow  rate  values  is  400  to  1000  in3/sec  which 
corresponds  to  a  recoil  velocity  range  of  70  to  170  in/sec.  The  Cd  range  for  these  conditions  is  0.70  to 
0.90.  The  trajectories  are  monotonic  and  show  signs  of  approaching  asymptotic  values  without  the 
presence  of  the  ‘knee’  as  shown  for  Zone  1  results.  As  before,  the  trajectory  for  the  largest  are  ratio  is  far 
removed  form  the  rest. 

The  upper  plot  on  Figure  26  contains  results  for  Zone  3  flow  rates.  The  range  of  flow  rate  values 
is  1000  to  2600  in3/sec  which  corresponds  to  a  recoil  velocity  range  of  170  to  420  in/sec.  The  Cd  range  for 
these  conditions  is  0.78  to  0.94.  Response  trajectories  are  as  before,  increasing  with  increasing  flow  rate, 
with  a  downward  shift  as  the  area  ratio  is  increased.  The  lower  plot  on  Figure  26  contains  results  for  Zone 
4  flow  rates.  The  range  of  flow  rate  values  is  2600  to  4400  in3/sec  which  corresponds  to  a  recoil  velocity 
range  of  420  to  700  in/sec.  The  Cd  range  for  these  conditions  is  0.85  to  0.96.  Comments  regarding  the 
responses  in  this  zone  are  the  same  as  above. 

To  sum  up  the  entire  CFD  analysis,  we  are  presenting  all  of  the  results  on  a  single  plot.  This  is 
presented  in  Figure  27.  The  flow  rate  range  is  0  to  4400  in3/sec  whereas  the  Cd  range  is  0.50  to  0.95.  For 
high  speed  flow  rates  (i.e.>2000  in3/sec)  the  Cd  range  is  between  0.85  and  0.95.  For  the  slower  flow  rates 
(500  to  2000  in3/sec)  the  Cd  range  is  between  0.75  and  0.92.  For  slow  flow  rates  (<  500  in3/sec)  the  Cd 
range  is  between  0.50  and  0.87.  When  the  trajectories  are  viewed  in  this  expanded  range  it  is  quite  clear 
that  a  strong  relationship  between  Cd  and  flow  rate  exists  containing  a  discernible  ‘knee’  around  the  500 
in3/secflow  rate.  Since  the  range  of  flow  rate  values  exercised  in  this  model  envelope  the  entire  range  of 
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Figure  27.  Cd  Results  for  all  Calculations 
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expected  values  and  enough  structure  has  been  demonstrated  for  each  area  ratio,  a  linear  interpolation 
of  Cd  between  flow  rate  values  is  warranted. 


Although  the  dependence  upon  area  ratio  is  not  as  severe,  it  is  quite  clear  that  as  this  value 
increases,  the  entire  range  of  Cd  values  decreases.  Unfortunately,  we  could  not  conduct  CFD  simulations 
for  the  very  large  area  ratios  since  the  time  and  cost  burden  would  be  prohibitive.  The  table  below 
presents  a  breakout  of  the  costs  for  each  area  ratio  exercised  in  these  simulations.  The  first  column 
contains  the  model  number  running  consecutively  from  1  to  19.  The  second  column  contains  the  radial 
values  of  the  control  rod  profile  whereas  the  third  contains  the  area  ratio  for  each  of  these  radial  values. 
The  fourth  column  contains  the  CPU  time  in  hours  for  the  33  simulations  (flow  rate  is  the  independent 
variable)  for  these  area  ratios  with  the  total  number  of  iterations  is  in  the  fifth  column.  The  last  column 
contains  the  time  per  iteration  in  seconds. 


GEOMETRIC 

MODEL 

NUMBER 

CONTROL 

ROD 

RADIUS 

(in) 

AREA 

RATIO 

CPU  TIME 
(hrs) 

NUMBER  OF 
ITERATIONS 

ITERATION 
RATE 
(secs  /  iter) 

1 

0.375 

11.08 

4.02 

14010 

2 

0.385 

11.57 

6.38 

15575 

1.47 

3 

0.395 

12.12 

5.81 

14740 

1.42 

4 

0.405 

12.75 

6.69 

16872 

1.43 

5 

0.415 

13.46 

4.39 

16657 

0.95 

6 

0.425 

14.27 

5.55 

15555 

1.28 

7 

0.435 

15.22 

8.02 

16519 

1.75 

8 

0.445 

16.32 

7.88 

16101 

1.76 

9 

0.455 

17.63 

6.73 

14074 

1.72 

10 

0.465 

19.21 

9.48 

15776 

2.16 

11 

0.475 

21.14 

8.61 

15933 

1.95 

12 

0.485 

23.55 

6.40 

14468 

1.59 

13 

0.495 

26.66 

8.54 

15087 

2.04 

14 

0.505 

30.82 

13.71 

18957 

2.60 

15 

0.515 

36.64 

18.24 

19136 

3.43 

16 

0.525 

45.38 

12.83 

18304 

2.52 

17 

0.535 

59.95 

35.95 

23023 

5.62 

18 

0.545 

89.12 

38.26 

15466 

8.91 

19 

0.555 

176.65 

277.53 

23610 

42.32 

Table  1 .  Summary  of  Cost  Burden  for  CFD  Analysis 
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As  we  examine  these  results  it  is  quite  clear  that  the  larger  area  ratios  are  quite  burdensome  in 
regard  to  execution  time.  In  addition  the  trend  is  non-linear  in  that  the  time  increases  by  a  factor  of  70 
from  the  smallest  to  the  largest  ratio  and  the  iteration  rate  increases  by  40.  The  final  ratio  reported  is 
176.65  which  correspond  to  a  control  rod  radius  of  0.555  in.  The  CPU  time  required  for  these  simulations 
was  over  1 1  days!  The  greatest  radius  of  the  current  control  rod  is  0.5625  which  occurs  at  the  very  end  of 
recoil  stroke.  The  area  ratio  for  this  rod  would  be  -700  and  the  number  of  elements  and  the  time  needed 
to  analyze  this  configuration  would  be  unfathomable.  So,  we  decided  to  use  an  extrapolation  method  to 
approximate  the  Cd  results  out  to  the  area  ratio  of  700. 

Figure  28  presents  the  calculations  (scatter  plots)  for  Cd  as  a  function  of  area  ratio  for  6  of  the  33 
flow  rate  values  used  in  the  CFD  analysis.  Color  coding  is  used  to  distinguish  between  different  flow 
rates.  In  the  upper  plot  the  results  are  presented  on  Cartesian  coordinates  whereas  in  the  lower  plot  a 
semi-log  graph  is  employed  with  the  common  logarithm  of  the  area  ratio  as  the  independent  variable.  The 
range  of  flow  rates  presented  bound  a  recoil  speed  of  5  to  600  in/sec.  The  area  ratio  axis  range  is  0  to 
200  for  the  upper  plot  and  10  to  700  for  the  lower  plot.  The  trend  as  a  function  of  area  ratio  is  similar  for 
all  flow  rates  with  a  decreasing  dependence  as  the  flow  rate  increases.  For  example,  in  the  upper  plot  the 
results  for  a  flow  rate  of  31  in3/sec  are  the  lowest  set  of  ‘data’  points  which  range  from  0.75  to  0.49.  The 
upper  set  of  ‘data’  ranging  from  0.95  to  0.88  is  for  flow  rate  of  3726  in3/sec.  The  dependence  upon  area 
ratio  becomes  less  severe  as  the  flow  rate  ratio  increases.  In  the  upper  plot  the  indicated  trend  is 
exponential  as  a  function  of  area  ratio.  When  the  common  logarithm  of  the  area  ratio  is  used  as  the 
independent  variable  the  points  tend  to  become  linear.  This  is  indicated  in  the  lower  plot.  The  lines  drawn 
within  each  set  of  ‘data’  merely  indicate  that  an  approximate  linear  trend  could  be  employed  to 
extrapolate  the  results  out  to  the  largest  area  ratio  value.  As  is  indicated  by  these  lines,  the  Cd  values  at 
an  area  ratio  of  700  could  be  reduced  by  as  much  as  0.15  from  the  lowest  calculated  value  at  area  ratio 
of  176.  This  response  characteristic  will  be  exploited  when  the  interpolation  algorithm  is  developed  for  the 
recoil  analysis  code. 

Brief  Explanation  of  the  Current  Recoil  Analysis  Code 

As  mentioned  earlier,  the  current  Benet  Labs  Recoil  Analysis  Code  employs  a  fixed  value  for  the 
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Figure  28.  Cd  Results  for  Various  Area  Ratios 
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discharge  coefficient  across  the  main  orifice.  Its  value  has  been  based  upon  those  published  in  technical 
literature  as  well  as  empirical  results  based  upon  firing  data.  The  most  applicable  values  are  usually 
between  0.90  and  0.95.  In  order  to  explain  the  implementation  of  the  variable  Cd  feature,  a  brief  ‘tour’  of 
the  code’s  structure  and  logic  flow  is  warranted.  The  initial  release  of  this  code  occurred  in  1981.  Since 
then,  it  has  been  modified  to  include  additional  features  and  more  elaborate  modeling  techniques.  A 
complete  rewrite  of  the  code  was  completed  in  the  early  1990’s,  after  which  several  other  modifications 
were  made  and  features  added. 

Figure  29  contains  a  skeletal  and  simplified  view  of  the  overall  code  via  standard  flowchart 
terminology.  The  INPUT  DATA  block  reads  overall  system  data  including  number  of  brakes  and 
recuperators,  weight  of  recoiling  parts,  bore  area,  friction  coefficients,  muzzle  brake  characteristics, 
integration  and  output  time  steps  as  well  as  total  analysis  time.  A  tabular  file  of  bore  pressure  versus  time 
is  input  as  well.  An  initial  call  is  made  to  the  RECUPERATOR  and  RECOIL  BRAKE  subroutines  to  read 
data  for  these  components.  Up  to  6  input  files  are  available  for  each.  The  recoil  brake  data  includes 
diameter  specifications  for  the  brake  cylinder  and  throttling  sleeve,  orifice  diameter,  temperature  and  bulk 
modulus  of  the  fluid,  several  discharge  coefficients  for  both  recoil  and  counter  recoil,  as  well  as  ramp  up 
distribution  features  for  partially  filled  fluid  chambers.  In  separate  tabular  files,  the  diameter  profile  of  the 
control  rod  and  the  characteristics  of  the  backflow  restrictive  areas  are  input  as  function  of  the  axial 
location  of  the  recoiling  parts.  A  file  containing  friction  factor  data  as  a  function  of  Reynolds  Number  (i.e. 
Moody  Chart)  is  input  and  is  used  in  calculating  the  pressure  drop  along  the  flow  length  to  and  from  the 
buffer  chamber.  Additionally,  files  for  the  density  and  viscosity  of  the  fluid  are  input  as  a  function  of  fluid 
temperature.  The  temperature  value  in  the  brake  input  file  will  determine  the  values  of  density  and 
viscosity  to  be  used  during  the  analysis.  The  recuperator  data  includes  relevant  diameter  and  length 
specifications  for  the  cylinder  and  rod,  preload  pressure  and  ratio  of  specific  heats  for  the  working  gas. 
Subsequently,  several  initial  values  within  the  SET  INITIAL  VALUES  block  are  set.  An  example  would  be 
the  values  of  the  integrator  coefficients  for  the  Adams-Bashfort  multistep  method  and  the  initial  values  for 
recoil  velocity  and  location. 
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Figure  29.  Overview  Flowchart  of  Benet  Labs  Recoil  Analysis  Code 
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The  dynamic  analysis  commences  in  the  START  INTEGRATION  LOOP.  The  ballistic  driving 
force  is  determined  based  upon  the  current  time  in  the  analysis  and  separate  calls  to  the  recuperator  and 
brake  subroutines  are  made  for  each  component  yielding  values  for  the  brake  and  recuperator  forces  at 
the  current  time  step.  Additional  forces  are  then  determined  (e.g.  friction,  etc)  from  which  the  recoil 
acceleration  for  the  current  time  step  is  calculated.  The  back  derivatives  are  updated  in  the  UPDATE 
BACK  DERIVATIVES  block  and  a  call  is  made  to  the  ODE  integrator  (CALL  ODE  INTEGATOR),  yielding 
the  current  values  for  recoil  velocity  and  travel.  The  PRINT  TIME  decision  test  yields  transient  values  for 
output  variables  when  applicable.  The  END  TIME  decision  test  either  returns  control  to  the  integration 
loop  for  continuation  of  another  time  step  or  ends  the  analysis.  Upon  analysis  termination,  maximum 
values  for  selected  output  variables  are  reported  to  a  maximum  values  file. 

Figure  30  contains  a  simplified  view  of  the  flow  within  the  recoil  brake  subroutine.  There  would  be 
multiple  calls  to  this  routine  within  each  time  step  dependent  upon  the  number  of  brakes  used  in  the 
analysis.  Upon  initial  entry  the  parameter  values  for  brake  ‘n’  are  placed  in  local  variables.  These  values 
include  brake  specification  dimensions,  current  flow  rates,  etc.  The  RECOIL  OR  C’RECOIL  decision 
diamond  determines  the  direction  of  travel.  There  are  separate  calculation  equations  for  each  direction. 
Dependent  upon  the  recoil  direction,  forces  and  pressures  are  determined  via  the  two  CALCULATE 
blocks.  Fluid  flow  rate  states  are  reset  for  the  current  brake  and  the  subroutine  returns  control  to  the  main 
program. 

Development  of  Variable  Discharge  Coefficient  Capability  in  the  Recoil  Analysis  Code 

Figure  31  contains  details  of  the  logic  within  the  branch  flow  portion  of  the  recoil  phase  calculation 
block.  As  mentioned  earlier,  the  flow  within  the  brake  during  recoil  follows  two  paths.  During  any  time  step 
the  value  of  the  outflow  rate  (in3/sec)  from  the  main  chamber  is  known  and  is  based  upon  the 
multiplication  of  recoil  speed  and  cross  section  area  of  the  main  fluid  chamber  of  the  brake  cylinder.  A 
portion  of  the  fluid  flows  through  the  annular  restrictor  defined  by  the  control  rod  and  the  main  orifice  while 
the  rest  flows  back  over  the  buffer  plug  into  the  buffer  chamber.  However,  a  closed  form  solution  to  the 
equations  yielding  pressure  for  these  flows  is  not  possible  since  one  branch  requires  the  use  of  friction 
factors  derived  from  a  Moody  Chart.  So  one  must  resort  to  estimating  the  flow  rate  split  (automatically 
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Figure  30.  Overview  Flowchart  of  Recoil  Brake  Subroutine 
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Figure  31.  Detail  Flowchart  of  Recoil  Pressure/Force  Calculations 
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done  within  the  routine)  and  calculate  the  resulting  pressures  based  upon  these  individual  flow  rates. 

Thus  the  initial  setting  of  an  expected  flow  rate  split  within  the  SET  INITIAL  SPLIT  FLOW  RATES  block  is 
accomplished.  The  logic  portion  USE  OF  VARIABLE  Cd  has  been  added  to  the  latest  version  of  the  recoil 
model.  Previous  versions  used  only  a  fixed  Cd  which  was  input  through  the  brake  data  file.  In  the  current 
model,  one  has  the  option  of  using  either.  When  the  variable  Cd  option  is  exercised,  the  VARIABLE  Cd 
SUBROUTINE  is  called  and  returns  the  current  value  for  Cd  based  upon  expected  flow  rate  through  the 
orifice  and  the  area  ratio,  which  is  a  function  of  the  location  of  the  recoiling  parts.  It  is  essentially  a  two 
parameter  lookup  table  using  linear  interpolation  for  the  flow  rate  and  linear  interpolation  of  the  common 
logarithm  of  the  area  ratio  to  calculate  values  between  data  points.  Upstream  fluid  pressure  is  calculated 
using  each  flow  rate  and  the  appropriate  equations  for  the  flow  along  that  path  yielding  two  values  for  the 
brake  pressure.  However,  the  physics  of  the  problem  dictates  that  both  pressures  must  be  the  same. 

Thus  the  two  are  compared  within  the  Po  ~  Pb  decision  diamond.  If  both  are  approximately  equal  (a  small 
residual  is  used  for  control)  then  the  flow  rates  along  each  path  are  correct  and  brake  force  is  calculated. 
If  not,  then  the  FLOW  SPLIT  ALGORITHM  block  is  accessed  and  based  upon  the  difference  in  pressure, 
the  flow  rate  along  each  path  the  is  adjusted  and  the  pressure  calculations  are  run  again.  This  loop  is 
traversed  until  the  difference  between  the  two  pressures  is  within  the  residual  value,  after  which  the  brake 
forces  are  determined  and  the  flow  volumes  within  both  the  front  and  buffer  chambers  are  updated. 

Figure  32  contains  the  flowchart  for  the  variable  Cd  subroutine.  The  data  passed  to  this  routine  is 
the  expected  flow  rate  and  the  area  ratio  at  the  current  location  of  the  recoiling  parts.  Upon  entry  a  test  is 
made  to  determine  if  this  is  the  first  pass  through  the  routine  (INITIAL  ENTRY  decision  diamond).  If  true, 
the  empirical  data  is  input  from  the  file  relating  Cd  to  flow  rate  and  area  ratio.  Subsequently,  the  entered 
values  (Q-in  &  A-in)  are  located  within  bounding  pairs  of  the  independent  variables.  From  here,  a  simple 
linear  interpolation  is  conducted  to  determine  the  Cd  value  to  be  used.  Control  is  returned  to  the  recoil 
brake  subroutine  with  the  applicable  value  for  Cd.  A  listing  of  the  FORTRAN  code  used  to  determine  the 
Cd  value  may  be  found  in  the  appendix. 


52 


Figure  32.  Detail  Flowchart  of  Variable  Cd  Subroutine 
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Results  Applied  to  Current  Data  from  XM360  Firing  Tests 

Comparison  of  Simulations  Using  Fixed  and  Variable  Discharge  Coefficients 

As  an  initial  estimate  of  the  effect  of  the  variable  coefficient  feature,  we  shall  conduct  simulations 
using  a  fixed  and  variable  Cd  and  compare  the  brake  pressure  results.  The  data  generated  in  the  CFD 
analysis  at  a  fluid  temperature  of  75°F  will  be  used  as  Cd  input.  The  ballistic  pressure  data  from  three 
120mm  rounds  will  be  used;  namely  M829A3,  M831A1  and  M865.  The  control  rod  is  the  as-designed 
version  for  a  recoil  travel  of  23.5  inches,  which  is  being  developed  for  the  Future  Combat  System  (FCS). 
Results  of  this  analysis  will  be  reported  in  pressure  versus  time  and  travel  plots  as  well  as  Cd  values 
plotted  against  time  and  travel. 

Brake  pressures  for  the  M829A3  round  are  shown  on  Figure  33.  The  upper  plot  shows  pressure 
versus  time  whereas  the  lower  plot  presents  pressure  versus  travel.  The  blue  line  is  the  response  for  a 
fixed  value  of  Cd ,  whereas  the  red  line  is  the  response  for  variable  Cd  values.  Peak  pressure  values  for 
both  are  about  4200  to  4500  psi.  During  the  initial  state  of  recoil  (time  <  15ms;  travel  <  4  inches)  both 
responses  track  each  other.  From  this  point,  the  variable  Cd  model  diverges  quickly  and  overshoots  the 
fixed  model  by  about  300  psi.  This  occurs  just  short  of  20  ms  at  a  travel  value  of  5  inches.  This  response 
separation  continues  up  to  about  60  ms  and  17.5  inches  of  travel.  At  this  point,  the  pressure  responses 
cross  each  other  and  the  variable  model  predicts  pressure  values  that  are  up  to  500  psi  lower  than  the 
fixed  model.  Total  recoil  stroke  for  the  variable  model  is  about  0.5  inches  shorter  than  the  same  for  the 
fixed  model  whereas  the  timing  has  only  a  3  ms  difference.  Cd  values  are  plotted  against  time  and  travel 
for  the  M829A3  round  on  Figure  34.  The  fixed  value  of  0.95  is  shown  along  with  the  values  determined 
within  the  dynamics  analysis  loop  in  the  recoil  analysis  routine.  The  variable  Cd  response  starts  quite  low 
but  then  rises  quickly  to  its  peak  value  of  0.94  within  10  ms  of  shot  initiation  and  2  inches  of  travel.  At  its 
peak,  the  Cd  variation  between  fixed  and  variable  is  only  0.01.  The  Cd  reduces  slightly  at  15  ms  (5  inches 
of  travel)  and  remains  at  this  level  until  about  30  ms  (10  inches  of  travel).  From  this  time  and  travel 
location,  the  Cd  reduces  very  slowly  to  about  0.86  at  60  ms  and  19  inches  of  travel.  It  then  reduces 
rapidly  eventually  to  a  value  of  0.45  at  1 1 0  ms  which  is  not  shown  on  the  plot. 
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Figure  33.  Simulated  Brake  Pressure  Comparison  for  M829A3  Round 
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Figure  34.  Simulated  Discharge  Coefficient  Comparison  for  M829A3  Round 
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Figure  35.  Simulated  Brake  Pressure  Comparison  for  M831A1  Round 
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Brake  pressures  for  the  M831 A1  round  are  shown  on  Figure  35.  The  upper  plot  shows  pressure 
versus  time  whereas  the  lower  plot  presents  pressure  versus  travel.  The  blue  line  is  the  response  for  a 
fixed  value  of  Cd ,  whereas  the  red  line  is  the  response  for  variable  Cd  values.  Peak  pressure  values  for 
both  are  about  3550  to  3850  psi.  During  the  initial  state  of  recoil  (time  <  15ms;  travel  <  3  inches),  both 
responses  track  each  other.  From  this  point,  the  variable  Cd  model  diverges  quickly  and  overshoots  the 
fixed  model  by  about  300  psi.  This  occurs  just  short  of  18  ms  at  a  travel  value  of  5  inches.  This  response 
separation  continues  up  to  about  60  ms  and  17.0  inches  of  travel.  At  this  point,  the  pressure  responses 
cross  each  other  and  the  variable  model  predicts  pressure  values  that  are  up  to  400  psi  lower  than  the 
fixed  model.  Total  recoil  stroke  for  the  variable  model  is  about  0.6  inches  shorter  than  the  same  for  the 
fixed  model,  whereas  the  timing  difference  is  negligible.  Cd  values  are  plotted  against  time  and  travel  for 
the  M831 A1  round  on  Figure  36.  The  fixed  value  of  0.95  is  shown  along  with  the  values  determined  within 
the  dynamics  analysis  loop  in  the  recoil  analysis  routine.  The  variable  Cd  response  starts  quite  low  but 
then  rises  quickly  to  its  peak  value  of  0.94  within  1 0  ms  of  shot  initiation  and  2  inches  of  travel.  At  its 
peak,  the  Cd  variation  between  fixed  and  variable  is  only  0.01.  The  Cd  reduces  slightly  at  17  ms  (4  inches 
of  travel)  and  remains  at  this  level  until  about  30  ms  (10  inches  of  travel).  From  this  time  and  travel 
location,  the  Cd  reduces  very  slowly  to  about  0.85  at  55  ms  and  1 6  inches  of  travel.  It  then  reduces  rapidly 
eventually  to  a  value  of  0.45  at  120  ms  which  is  not  shown  on  the  plot. 

Brake  pressures  for  the  M865  round  are  shown  on  Figure  37.  The  upper  plot  shows  pressure 
versus  time  whereas  the  lower  plot  presents  pressure  versus  travel.  The  blue  line  is  the  response  for  a 
fixed  value  of  Cd,  whereas  the  red  line  is  the  response  for  variable  Cd  values.  Peak  pressure  values  for 
both  are  about  2150  to  2350  psi.  During  the  initial  state  of  recoil  (time  <  12ms;  travel  <  3  inches),  both 
responses  track  each  other.  From  this  point,  the  variable  Cd  model  diverges  quickly  and  overshoots  the 
fixed  model  by  about  200  psi.  This  occurs  just  short  of  19  ms  at  a  travel  value  of  5  inches.  This  response 
separation  continues  up  to  about  75  ms  and  17.0  inches  of  travel.  At  this  point,  the  pressure  responses 
cross  each  other  and  the  variable  model  predicts  pressure  values  that  are  up  to  200  psi  lower  than  the 
fixed  model.  Total  recoil  stroke  for  the  variable  model  is  about  0.1  inches  shorter  than  the  same  for  the 
fixed  model  whereas  the  timing  difference  is  1 .5  ms.  Cd  values  are  plotted  against  time  and  travel  for  the 
M831 A1  round  on  Figure  38.  The  fixed  value  of  0.95  is  shown  along  with  the  values  selected  within  the 
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Figure  36.  Simulated  Discharge  Coefficient  Comparison  for  M831A1  Round 
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dynamics  analysis  loop  in  the  recoil  analysis  routine.  The  variable  Cd  response  starts  quite  low  but  then 
rises  quickly  to  its  peak  value  of  0.93  within  8  ms  of  shot  initiation  and  2  inches  of  travel.  At  its  peak,  the 
Cd  variation  between  fixed  and  variable  is  only  0.02.  The  Cd  reduces  slightly  to  about  0.90  at  17  ms  (4 
inches  of  travel)  and  remains  at  this  level  until  about  38  ms  (1 0  inches  of  travel).  From  this  time  and  travel 
location,  the  Cd  reduces  very  slowly  to  about  0.83  at  65  ms  and  17  inches  of  travel.  It  then  reduces  rapidly 
eventually  to  a  value  of  0.45  at  140  ms  which  is  not  shown  on  the  plot. 

In  general,  the  comparison  between  the  fixed  and  variable  models  was  consistent  for  all  three 
round  types  considered.  The  variable  model  produced  peak  pressure  values  slightly  greater  than  the 
fixed.  The  pressure  response  for  the  variable  model  exceeded  the  same  for  the  fixed  for  about  one-half  of 
the  recoil  cycle  in  time  and  three-quarters  in  travel.  Subsequently,  the  variable  model’s  response  was  less 
than  the  fixed.  The  maximum  Cd  value  for  the  variable  model  was  slightly  less  than  the  expected  value  of 
0.95  and  reduced  in  a  similar  manner  as  time  and  recoil  stroke  progressed.  The  final  values  for  Cd  always 
became  quite  small  as  recoil  speed  decreased  and  area  ratio  increased. 

Comparison  of  Test  Data  to  Predictions  from  Simulations 

In  order  to  determine  the  benefit  in  using  this  new  analysis  algorithm,  we  shall  compare  results 
against  brake  pressure  data  collected  during  firing.  During  the  various  ATD  tests  of  the  XM360  weapon 
there  were  occasions  during  which  several  shots  were  fired  in  a  relatively  short  period  of  time.  For 
example,  the  Target  Impact  Dispersion  (TID)  test  dictated  that  the  firing  cadence  be  as  quick  as  possible. 
As  a  norm,  these  tests  usually  required  the  firing  of  at  least  1 0  rounds  to  get  meaningful  data.  During  the 
ATD-1  tests  three  round  types  were  fired  in  this  manner  during  a  2  day  period.  We  shall  use  these  results 
to  assess  the  worth  of  the  variable  Cd  as  compared  to  a  fixed  Cd  analysis  model.  As  an  indication  of  the 
robustness  of  the  recoil  system,  we  shall  compare  brake  pressure  data  for  several  rounds  on  a  single  plot 
and  superimpose  the  average  response  of  the  rounds. 

The  first  of  these  is  shown  on  Figure  39  on  which  the  brake  pressure  data  response  for  9  shots  of 
the  M829A3  round.  The  upper  plot  contains  the  data  for  the  upper  left  brake  whereas  the  lower  plot  has 
the  data  for  the  lower  right  brake.  The  ‘dust  particle’  plot  (points  only)  is  the  data  results  for  all  9  shots 
whereas  the  heavy  line  plot  is  the  average  of  the  9  shots.  The  similarity  of  the  response  for  each  brake  is 
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Figure  39.  ATD-1  Test  Results  Brake  Pressure  M829A3 
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noteworthy.  Additionally,  the  variation  from  shot  to  shot  is  minimal  indicating  stability  in  the  design  of  the 
components.  The  meandering  of  the  ‘dust  particles’  from  the  average  is  only  indicated  during  the  early 
portions  of  the  response  and  just  beyond  the  peak  level.  Subsequent  to  this,  pressure  values  from  the 
individual  responses  are  indeterminable  from  the  average.  Figures  40  and  41  contain  the  same  type  of 
data  for  the  M831 A1  and  M865  rounds.  For  these  rounds,  the  within  round  variation  in  the  pressure 
response  is  slightly  larger,  especially  early  in  the  cycle.  However,  the  shot  to  shot  stability  is  quite  good. 
The  average  data  for  these  round  types  will  be  used  as  comparison  with  the  simulation  runs. 

These  test  firings  were  conducted  in  August  of  2004  at  Aberdeen  Test  Center  where  the  ambient 
temperature  range  is  usually  between  90°F  and  100°F.  Since  we  do  not  have  fluid  nor  brake  cylinder 
temperature  data  for  these  tests,  a  reasonable  estimate  would  be  between  100°F  to  125°F  and  possibly 
higher  depending  upon  the  intensity  and  duration  of  radiant  solar  heat.  One  would  think  that  this  would 
produce  a  marginal  change  in  results,  but  the  fact  is  the  viscosity  of  the  fluid  at  100°F  is  about  60  percent 
of  its  value  at  75°F.  The  density  change,  however,  is  negligible.  An  initial  set  of  simulations  and 
comparisons  were  completed  using  the  ‘data’  from  the  Cd  tables  developed  using  a  fluid  temperature  of 
75°F  and  the  results  were  not  bad  but  the  peak  pressure  was  over-predicted  by  200  to  700  psi  depending 
upon  round  type. 

At  this  point  in  the  project,  we  decided  to  reevaluate  the  Cd  calculations  at  an  elevated 
temperature  of  1 00°F  to  determine  its  sensitivity  to  temperature.  The  complete  list  of  CFD  models  were 
rerun  at  a  fluid  temperature  of  100°F  and  the  results  collected  for  comparison.  The  CFD  results  for  area 
ratios  of  23.5  and  12.7  are  shown  on  Figure  42.  On  this  figure  Cd  is  plotted  against  flow  rate  for  the  fluid 
temperature  values  of  75°F  and  100°F.  As  indicated  the  distribution  for  the  higher  fluid  temperature  model 
is  slightly  to  significantly  greater  than  the  lower  temperature.  For  the  smaller  area  ratio  value  (1 2.7)  which 
would  occur  early  in  the  cycle  at  high  flow  rate  values,  the  Cd  values  are  0.01  to  0.03  units  greater  at  the 
higher  temperature.  At  the  larger  area  ratio  (23.5)  the  Cd  values  are  greater  at  the  higher  temperatures 
but  eventually  converge  at  the  highest  flow  rates.  The  additional  benefit  gleaned  from  running  this 
additional  CFD  analysis  was  that  we  could  extrapolate  to  other  proximate  temperatures  without  the  need 
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Figure  40.  ATD-1  Test  Results  Brake  Pressure  M831A1 
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Figure  41.  ATD-1  Test  Results  Brake  Pressure  M865 
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Figure  42.  CFD  Calculations  @  Various  Fluid  Temperatures 
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to  run  the  complete  CFD  model  at  the  expected  temperature.  For  the  comparisons  detailed  below,  we 
extrapolated  the  Cd  values  to  125°F  using  the  CFD  ‘data’  from  the  two  temperature  dependent  runs  (75°F 
and  100°F)  and  assumed  a  linear  response  beyond  the  temperature  range.  We  cannot  justify  the  linear 
extrapolation  until  additional  analyses  are  conducted,  however,  for  this  exercise  this  should  be 
acceptable. 

Figures  43  through  45  contain  graphical  results  comparing  the  test  data  and  simulations.  The 
graphs  on  these  figures  present  brake  pressure  versus  time  in  the  upper  plot  and  pressure  versus  travel 
in  the  lower  plot.  Within  each,  three  graphs  entitled  TEST  DATA,  VARIABLE  Cd  @  125°F  and  FIXED  Cd 
@  125°F  are  indicated.  The  assumed  value  for  the  fluid  temperature  was  125°F,  thus  the  table  of  Cd 
values  had  to  be  extrapolated  to  125°F  using  the  values  at  75°F  and  100°F.  The  test  data  used  for 
comparison  is  the  average  response  for  the  series  of  shots  as  indicated  previously.  The  Cd  value  used  for 
the  fixed  model  was  0.95. 

Results  for  the  M829A3  round  are  shown  on  Figure  43.  As  indicated,  the  simulated  responses  for 
either  fixed  or  variable  Cd  are  quite  similar.  Before  peak  pressure  is  reached  (up  to  17ms  and  5  inches  of 
travel)  the  results  from  both  models  are  essentially  superimposed.  However,  at  peak  pressure  the 
variable  model  predicts  a  pressure  value  that  is  slightly  less  than  its  fixed  Cd  counterpart.  Additionally,  the 
“sharpness  factor”  which  is  the  plateau-type  flat  is  characteristic  for  both  models.  Early  in  the  cycle,  test 
data  indicates  a  response  that  is  greater  than  either  model.  However,  at  peak  pressure  the  data  response 
is  slightly  lower  than  either,  its  timing  somewhat  later  (20  ms  and  6  inches)  and  its  duration  much  broader 
without  the  “sharpness  factor”  that  is  indicated  in  the  simulations.  Beyond  peak  pressure  and  up  to  60  ms 
and  15  inches  of  travel,  both  models  over  predict  the  data  by  as  much  as  300  to  500  psi  depending  upon 
the  model.  The  response  from  the  variable  Cd  model  is  more  consistent  with  the  data  than  the  fixed  Cd 
model  during  this  portion  of  the  cycle.  At  90  ms,  both  simulation  models  and  the  data  converge  at  about 
800  psi  and  the  end  of  cycle  occurs  at  1 1 0  ms  for  both  models  and  the  data.  Maximum  travel  for  variable 
model  exceeds  the  data  by  about  one-half  inch,  whereas  for  the  fixed  model  the  difference  is  about  one 
inch. 
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Figure  43.  Comparison  of  Brake  Pressure  for  M829A3  Round 
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Figure  45.  Comparison  of  Brake  Pressure  for  M865  Round 
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Results  for  the  M831 A1  round  are  shown  on  Figure  44.  In  a  manner  that  is  similar  to  the  previous 
round,  the  simulated  responses  for  either  fixed  or  variable  Cd  are  quite  similar.  Both  track  exactly  prior  to 
the  achievement  of  the  peak  pressure  value.  The  “sharpness  factor”  exists  at  peak  pressure  with  a  value 
that  exceeds  the  data  by  about  400  psi.  The  prediction  from  the  variable  Cd  model  is  slightly  greater  than 
for  the  fixed.  During  certain  portions  in  the  cycle  beyond  peak  pressure,  the  prediction  from  variable 
model  exceeds  the  fixed  model  by  a  few  hundred  psi.  During  the  decay  portion  both  simulation  responses 
converge  at  100  ms.  The  total  cycle  time  for  both  models  and  the  data  is  around  120ms.  The  results 
plotted  against  travel  indicate  the  same  trends  as  reported  above.  Both  simulation  models  slightly  over 
predict  pressure  for  the  complete  length  of  travel.  In  addition,  the  predicted  maximum  travel  slightly 
exceeds  that  of  the  data  by  about  one-half  inch. 

Results  for  the  M865  round  are  shown  on  Figure  45.  With  the  exception  of  the  first  1 5  ms  of  data, 
the  similarities  and  discrepancies  indicated  in  the  previous  results  are  the  same  for  this  round.  Both 
simulation  models  over  predict  the  data  by  at  most  a  few  hundred  psi  with  the  differences  from  variable 
model  slightly  greater  than  the  fixed.  However,  maximum  travel  as  predicted  by  both  models  is  the  same 
as  the  data.  Regarding  the  early  portion  of  the  cycle  and  recalling  that  the  data  and  average  for  the  10 
M865  rounds  were  shown  on  Figure  42,  the  ‘dust’  plots  of  the  individual  shots  during  this  portion  of  the 
cycle  were  rather  ‘chaotic’.  This  resulted  in  an  average  that  captured  the  erratic  nature  of  the  data  as  well. 

On  Figure  46  the  calculated  Cd  values  (determined  from  the  flow  rates  in  recoil  analysis  code)  for 
all  three  round  types  as  well  as  the  value  for  the  fixed  Cd  model  are  presented.  The  upper  plot  contains 
these  values  potted  against  time  whereas  on  the  lower  plot  these  results  are  plotted  against  travel.  The 
range  of  values  for  the  variable  Cd  is  0.50  to  0.98.  As  indicated,  the  Cd  values  achieve  their  maximum  of 
0.98  early  in  the  recoil  cycle.  For  the  first  4  inches  of  travel  and  20  ms  in  time,  these  sustained  levels  are 
greater  than  the  previously  assumed  value  of  0.95.  Subsequently,  the  trend  decreases  from  their 
maximum  values  to  a  sustained  level  ranging  from  0.91  to  0.93  for  a  significant  portion  of  the  cycle.  An 
abrupt  ‘knee’  in  all  three  plots  is  shown  to  occur  towards  the  end  of  the  cycle  when  20  to  22  inches  of 
travel  has  been  achieved.  Overall,  the  Cd  distribution  contains  areas  of  distinct  oscillations.  This  is  most 
likely  due  to  the  nature  of  the  fitting  function  used  for  finding  the  Cd  values  between  independent  variable 
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Figure  46.  Comparison  of  Calculated  Cd  Values 
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points.  Recalling  Figure  28,  on  which  the  Cd  response  as  a  function  of  area  ratio  is  shown,  the  trend  is 
linear  on  the  semi-logarithmic  plot,  however,  the  individual  ‘data’  points  meander  above  and  below  the 
trend  line.  The  fitting  function  used  in  the  analysis  uses  a  linear  fit  between  area  ratio  points  thus  the  Cd 
value  would  tend  to  oscillate  above  and  below  the  trend  line.  The  remedy  would  be  to  use  the  trend  line  to 
predict  the  Cd  values  or  reduce  the  number  of  area  ratio  points  in  the  table  of  values. 

Discussion  of  the  Significance  of  the  Fixed  and  Variable  Discharge  Coefficient  Method 

Although  the  results  did  not  indicate  predictive  capabilities  that  are  much  better  than  we  currently 
have  using  an  assumed  and  fixed  Cd,  we  did  demonstrate  that  the  using  off-line  CFD  results  as  ‘data’  to 
the  recoil  analysis  simulations  does  lend  itself  to  results  that  are  comparable  to  our  current  capabilities. 
Since  we  have  been  conducting  recoil  design  and  analysis  for  at  least  25  years  at  Benet  Labs,  we 
possess  the  experience  needed  to  estimate  appropriate  Cd  values  from  the  myriad  of  test  data  and 
simulation  runs  generated  during  this  time.  Since  the  characteristics  of  the  flow  path  from  the  annular 
chamber  of  the  brake  and  through  the  orifice  is  similar  to  all  previous  models  that  have  been  built  and 
tested,  the  choice  of  a  Cd  value  of  0.95  for  the  XM360  system  is  consistent  and  renders  very  accurate 
predictions  of  test  results.  In  addition,  when  the  Cd  values  from  CFD  analysis  and  dynamic  simulations 
are  compared  to  the  extrapolated  values  from  firing  data  (Figures  16-18),  the  distribution  is  similar.  The 
Cd  value  is  large  during  early  phases  of  recoil  (0  to  5  in),  and  then  decreases  slowly  as  recoil  stroke 
approaches  its  maximum  value. 

However,  there  are  still  other  considerations.  For  example,  the  accuracy  of  the  test  data  is 
questionable,  especially  during  the  early  portion  of  the  cycle  when  speeds,  flow  rates  and  area  ratios  are 
changing  rapidly.  The  data  from  Figures  43  through  45  indicated  a  rather  benign  peak  in  the  pressure 
response  whereas  the  simulations  show  a  sharp  and  distinct  peak.  Data  was  extracted  from  a  gage 
located  in  the  brake  cylinder’s  the  fill  port  which  is  not  embedded  directly  in  the  fluid  and  at  a  90°  offset 
from  the  direction  of  flow,  thus  the  recorded  data  may  not  be  commensurate  with  the  average  pressure  in 
the  stream.  Another  possible  contributor  is  the  elastic  response  of  the  cylinder  structure  that  may  have 
affected  the  accuracy  of  the  recorded  pressure.  Currently,  this  is  our  best  method  for  determining 
pressure  response  and  may  continue  to  be  for  some  time. 
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The  real  worth  of  this  exercise  may  be  in  determining  the  relationship  between  pressure  and  flow 
rate  for  flow  paths  and  orifice  geometries  that  are  not  consistent  with  our  current  models.  For  example,  in 
counter  recoil,  the  main  orifice  contributes  to  the  generation  of  counter  recoil  forces  in  a  manner  that  is 
similar  to  recoil,  however,  the  orifice  shape  that  is  presented  to  the  reverse  flow  is  ‘flange-like’  with  a 
sharp  taper  to  the  throat.  This  feature  is  much  different  than  in  the  recoil  direction.  Again  we  have  been 
using  Cd  values  derived  from  firing  data  (albeit,  we  could  not  get  reasonable  force  nor  kinematics  data 
until  the  XM360  program)  and  have  found  the  value  to  be  somewhat  lower  than  it  is  in  the  recoil  direction. 
However,  by  using  the  same  CFD  models  developed  for  the  recoil  direction,  the  Cd  values  in  counter 
recoil  may  be  determined  by  reversing  the  inlet  and  outlet  boundaries  and  the  direction  and  magnitude  of 
flow  rate.  The  Bernoulli  Equation  and  the  calculation  of  Cd  values  in  spirit  is  the  same  in  counter  recoil  as 
it  is  in  recoil.  Also,  the  flow  from  the  buffer  chamber  in  counter  recoil  is  considerably  more  complicated 
than  in  recoil,  moreover,  it  has  been  shown  through  unsteady  CFD  analysis  that  cavitation  bubbles  that 
form  during  recoil  collapse,  and  in  some  cases  pass  over  the  buffer  plug  and  into  the  rear  cavity  [8].  So,  in 
a  manner  similar  to  the  methods  presented  in  this  work,  we  may  use  CFD  analysis  to  predict  the 
relationship  between  flow  rate  and  fluid  pressure  in  the  counter  recoil  direction  for  flows  emanating  from 
the  buffer  chamber. 

Discussions  and  Conclusions 

Summary  of  Analysis  Philosophy  and  Results 

As  detailed  in  the  introduction,  the  motivation  driving  the  conduction  of  this  study  is  the  need  to 
more  accurately  model  the  flow  variables  within  the  chambers  of  recoil  brakes  in  a  dynamic  environment. 
Since  the  recoil  brake  provides  greater  than  ninety  percent  of  the  retarding  force  in  recoil,  it  is  essential 
that  the  most  accurate  method  of  modeling  this  component  is  employed.  Prior  to  this  work,  several 
assumptions  regarding  the  efficiency  of  the  brake’s  main  orifice  were  taken  from  published  values  for 
restrictors  of  similar  geometry.  Neither  the  fluid  characteristics  nor  the  flow  regime  was  factored  into  these 
values.  In  addition,  a  method  of  ‘backing  out’  these  values  from  field  generated  test  data  proved  to  be 
somewhat  subjective  as  well.  We  concluded  that  a  better  method  could  be  devised  by  combining  the 
results  of  an  off-line  and  more  detailed  CFD  analysis  to  the  simplified  recoil  brake  model  currently  used 
for  analysis.  There  have  been  attempts  at  using  a  full-blown  CFD  dynamic  analysis  model  of  recoil  brake 
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performance  but  due  to  the  3D  nature  of  the  flow  and  the  need  for  fine  grid  meshing,  a  real-time  analysis 
is  quite  difficult  to  achieve.  However,  the  use  of  parallel  processing  and  higher  speed  super  computers 
should  dramatically  speed  up  analysis  turn-around  times.  We  devised  a  method  of  modeling  the  intricate 
flow  paths  within  a  recoil  brake  and  calculate  the  orifice  efficiency  using  a  CFD  flow  model,  the  results  of 
which  were  tabulated  against  flow  areas  and  rates.  Although  this  required  several  CPU  hours  and  weeks 
to  conduct  these  studies,  once  completed  the  results  were  used  repeatedly  as  ‘data’  for  the  one  degree  of 
freedom  Benet  Labs  recoil  model. 

Unfortunately,  the  viscosity  value  of  the  fluid  (being  a  strong  function  of  fluid  temperature)  greatly 
influenced  the  Cd  values  calculated  in  the  CFD  analysis.  A  new  table  of  values  had  to  be  generated  in 
order  to  best  model  the  firing  data  since  the  fluid  temperature  of  the  data  was  approximately  1 1 0  F  to 
125  F  and  not  the  75°F  used  in  the  initial  CFD  analysis.  One  may  postulate  that  a  better  independent 
variable  may  be  the  Reynolds  Number  of  flow  through  the  orifice  since  fluid  properties  (at  any 
temperature),  flow  rate  and  geometric  structure  are  intrinsic  to  this  value.  To  determine  the  worth  of  this 
hypothesis  we  conducted  a  cursory  CFD  analysis  for  a  single  area  ratio  (12.12)  with  fluid  at  various 
temperatures  and  the  same  list  of  flow  rates  expected  during  actual  firing  of  the  weapon.  Seven  values  for 
fluid  temperature  along  with  their  comparable  values  for  specific  gravity  and  dynamic  viscosity  were  used 
along  with  sixteen  flow  rates.  Thus,  112  separate  runs  comprise  this  mini  analysis.  The  results  are 
presented  on  Figure  47  in  the  form  scatter  plots  of  Cd  versus  Reynolds  Number  using  two  axes  types.  In 
the  upper  plot  the  Cd  values  are  plotted  on  a  linear  scale  whereas  the  Reynolds  Number  is  shown  on  a 
common  logarithm  scale.  In  the  lower  plot  both  Cd  and  Reynolds  Number  values  are  shown  on  a  common 
logarithm  scale.  The  structure  of  these  results  is  quite  striking  and  indicates  extreme  tractability  in  that  the 
results  tend  to  meld  into  each  other  as  the  analysis  transitions  from  one  temperature  to  another.  The  grey 
symbols  are  results  at  fluid  temperature  of  -50°F  for  which  the  Reynolds  number  range  is  ~5  to  -300.  The 
cyan  symbols  at  the  other  extreme  are  for  a  temperature  of  250°F  for  which  the  Reynolds  number  range 
is  -8500  to  -600,000.  Thus,  the  entire  range  for  Reynolds  Number  is  5  to  600,000.  The  results  as 
presented  on  the  Log-Log  plot  indicates  more  linearity  for  the  lower  values  of  Reynolds  Number  as  well  a 
flatter  response  at  upper  values.  All  in  all,  it  is  the  same  set  of  data  and  the  reason  for  presenting  on  both 
axes  is  to  show  that  a  rather  benign  fitting  function  could  be  devised  to  extrapolate  values  between  ‘data’ 
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points,  yielding  very  acceptable  results.  We  shall  continue  exploring  this  idea  at  different  area  ratios.  In 
fact,  we  may  discover  that  since  area  ratio  is  mildly  contained  in  the  Reynolds  Number  (i.e.  the  hydraulic 
radius  is  contained  in  the  numerator),  all  models  may  meld  together.  This  would  quite  serendipitous,  in 
that  a  single  ‘data’  file  of  Cd  values  with  only  Reynolds  Number  as  the  independent  variable  could  be 
devised.  Since  the  Benet  Labs  recoil  model  requests  fluid  temperature  the  Reynolds  Number  may  be 
calculated  within  the  dynamic  loop  and  the  appropriate  Cd  value  found.  The  burden  would  only  be  the 
time  needed  to  set  up  and  extract  results  from  the  bank  of  runs  (which  is  minimal)  and  the  CPU  time 
needed  to  conduct  these  runs  (which  for  some  geometries  is  quite  extensive). 

When  applied  to  current  data  from  XM360  gun  firings,  this  variable  Cd  method  proved  to  be  as 
good  a  predictor  of  the  brake  pressure  data  as  the  previous  fixed  Cd  model.  For  some  regimes  of  the  flow, 
the  variable  model  more  accurately  predicted  pressure  values  than  did  the  fixed  and  vice  versa  for  other 
regimes.  In  any  event,  both  models  predict  pressure  values  that  are  quite  close  to  reality,  however,  if  we 
use  the  variable  model  we  no  longer  need  to  ‘guess’  at  Cd  values  which  has  always  been  problematic  at 
best.  In  fact,  several  other  restrictor  geometries  for  which  Cd  data  is  unavailable  could  be  analyzed  prior 
to  the  incorporation  into  brake  design  or  simulation.  The  only  downside  is  the  need  for  detailed  CFD 
analysis  to  determine  Cd  values  prior  to  conducting  recoil  simulations. 

Model  Enhancement  for  Future  Endeavors 

Since  this  methodology  proved  to  be  useful  for  the  recoil  portion  of  the  cycle  it  would  seem 
lucrative  to  conduct  similar  studies  for  the  counter  recoil  portion.  The  main  flow  orifice  which  controls 
brake  pressure  in  recoil  does  the  same  in  counter  recoil.  However,  the  geometry  and  cross  section  area 
are  not  the  same  nor  is  the  range  of  flow  rates.  (In  counter  recoil  the  velocity  barely  reaches  100  in/sec 
whereas  in  recoil  the  maximum  speed  is  -700  in/sec.)  A  separate  table  of  values  for  Cd  would  have  to  be 
included  in  the  recoil  model  but  the  method  for  its  incorporation  and  use  is  the  same  as  it  is  for  recoil. 

In  addition,  there  are  other  restrictive  flow  paths  through  which  the  fluid  flows  in  counter  recoil. 
The  flow  from  the  buffer  pocket  over  and  through  the  buffer  plug  is  quite  complicated.  There  are  two 
parallel  paths  through  which  the  flow  is  apportioned.  Currently,  we  use  several  coefficients  to  model 
abrupt  entrance  and  exit  effects  as  well  as  a  dedicated  Moody  Diagram  to  determine  friction  factors  along 
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the  flow  length  between  the  external  profile  of  the  buffer  plug  and  the  circumferential  surface  of  the 
throttling  sleeve.  This  path  has  tapers  and  abrupt  changes  in  diameter  as  well  as  circumferential  grooves 
to  control  the  back  pressure  in  the  buffer  pocket.  In  fact,  it  is  our  experience  that  modeling  flow  in  these 
regions  is  quite  a  bit  more  difficult  than  in  the  recoil  direction.  The  various  entrance  and  exit  coefficients 
used  in  the  current  model  have  been  determined  based  upon  the  results  of  test  data  in  counter  recoil  and 
have  been  set  to  best  match  the  data.  Since  the  total  recoil  time  is  a  highly  sensitive  value  for  the  XM360 
program  and  the  fact  that  counter  recoil  time  comprises  the  bulk  of  cycle  time  it  would  be  quite  beneficial 
to  be  able  to  model  flow  in  this  regime  as  well  as  predicting  terminal  velocity  and  total  cycle  time  prior  to 
test  firing. 

Recommendations  for  Future  Work 

As  with  most  comprehensive  analyses,  the  initial  findings  usually  beget  the  question  of  how 
accurate  are  the  results  or  the  fidelity  of  the  models  to  reality.  Earlier  we  set  the  ground  rules  for  the  CFD 
analysis  to  be  axisymmetric,  incompressible  and  steady  state.  We  choose  these  characteristics  since  the 
flow  equations  during  any  time  step  within  the  dynamics  loop  of  the  recoil  model  are  those  for  quasi¬ 
steady  state  conditions.  In  addition,  these  model  characteristics  allowed  the  conduction  of  manageable 
and  cost  effective  analyses  from  which  timely  results  were  extracted.  As  a  follow-on,  it  may  be  beneficial 
to  reproduce  the  pressure  and  flow  rate  relationships  using  another  more  detailed  CFD  code  namely 
Fluent®.  In  addition,  the  non  included  characteristics  of  asymmetry,  compressibility,  and  unsteady  motion 
could  be  modeled  using  both  CFDesign®  and  Fluent®  to  determine  their  effect.  In  addition,  Fluent®  could 
be  used  to  model  multi-phase  cavitated  flow  and  more  complicated  fluid  properties  that  vary  as  a  function 
of  temperature. 

Since  there  are  flow  paths  in  the  counter  recoil  direction  which  are  more  complicated  than  those 
in  recoil,  the  application  of  the  CFD  analysis  for  these  paths  should  be  conducted  and  the  results  used  in 
the  recoil  analysis  code.  However,  the  results  may  take  on  a  different  flavor.  For  the  recoil  direction,  the 
CFD  ‘data’  was  applied  to  the  steady-state  flow  equation  to  determine  Cd  values  used  in  the  recoil 
analysis  code.  In  counter  recoil,  since  there  are  flow  lengths  as  well  as  tapered  restrictors,  it  may  be 


79 


easier  to  use  the  pressure  versus  flow  rate  results  directly  in  the  recoil  model.  Currently,  we  use  entrance 
coefficients  and  a  modified  Moody  Diagram  to  determine  orifice  efficiency  and  flow  length  friction  factors 
to  achieve  the  pressure  drop  along  any  restrictive  path.  A  direct  table  of  values  relating  flow  rate  and 
recoil  position  to  back  pressure  would  alleviate  the  need  to  use  the  friction  diagrams  and  artificial 
entrance  coefficients. 

As  recounted  earlier,  the  time  required  to  have  a  representative  set  of  meaningful  data  from  the 
CFD  models  was  quite  extensive.  This  was  due  to  many  factors.  First,  the  use  of  the  CFDesign  code 
required  an  initial  ‘learning  curve’  since  these  authors  had  no  experience  in  its  use.  Second,  a  priori 
knowledge  in  the  expected  results  dictated  a  conservative  approach  to  the  grid  meshing  of  the  flow 
regime  as  well  as  the  number  and  spacing  of  independent  variable  values.  Since  we  now  have 
considerable  knowledge  in  these  issues,  design  of  experiments  (DOE)  techniques  may  be  employed  for 
future  use  of  this  methodology.  DOE  can  be  used  in  a  sequential  approach,  starting  with  simple  models 
and  adding  (augmenting)  as  needed.  This  could  have  reduced  the  number  of  runs  required.  DOE  has 
built  in  capabilities  within  the  ANOVA  tool  sets  to  tell  you  if  you  need  more  data  to  get  a  better  fit. 

Currently,  the  design  of  a  brake  system  begins  with  the  use  of  a  simple  spreadsheet  code,  the 
input  of  which  comprises  operating  characteristics  along  with  other  fixed  and  limiting  values  for  geometric 
structure  and  fluid.  The  orifice  coefficient  is  one  of  the  fixed  values  required.  We  have  always  used  an 
appropriate  value  determined  from  experience.  Output  from  the  design  code  contains  the  remaining 
geometry  values  for  brake  and  recuperator  diameters  and  lengths  as  well  as  a  tabulation  of  the  control 
rod  diametrical  profile.  We  treat  this  table  of  values  as  a  starting  point  for  the  design.  The  analysis  code 
uses  these  values  to  ‘refine’  the  profile  for  optimization  across  several  operating  parameters  and 
dimensional  tolerances.  Prior  to  this  step  and  after  the  basic  design  of  the  brake  is  known,  the  detailed 
CFD  analysis  would  be  used  to  determine  expected  Cd  values.  In  the  very  least,  basic  research  in  this 
area  should  continue  irrespective  of  any  design  program.  Several  different  restrictor  characteristics  could 
be  modeled  to  determine  the  most  efficient  and  optimum  regarding  the  size  and  weight  tradeoff  which  is 
always  at  the  forefront  of  any  weapon  design.  Additionally,  our  modeling  capabilities  in  recoil  would  rise  to 
the  next  level. 
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APPENDIX 


Matlab®  Cd  Calculator  Routine 


o, 

o 

%....  Initial  code  to  Determine  Cd  from  data  (1/30/07) 

%....  This  version  computes  Reynolds  # 

%....  load  1:  Input  test  data  from  firing  test 
%....  load  2:  Brake  geometric  data  from  drawing 

%....  load  3:  Travel  /  Velocity  'DATA'  from  Sigma  Plot  Analysis 

o, 

o 

%....  Contains  0/P  for  SP  plotting 

O, 

o 

load  Rnd_113.dat; 
load  BL_46569.dat; 

o, 

o 

%....  Set  poly  order  for  poly  fits  for  data 

o, 

o 

Nx  =  3; 

Np  =  20; 

Q, 

O 

%....  O/P  file  for  SP  plotting 

O, 

o 

filenamel  =  ' Rnd  113_L  R.dat'; 

O, 

o 

o, 

o 

%....  set  data  into  local  vectors 

o, 

o 

t  =  Rnd_l 1 3 ( : , 1 ) ; 
x  =  Rnd  113  (  :  , 2)  ; 
pi  =  Rnd_113 ( : , 3) ; 

pr  =  Rnd_113 ( : , 4) ; 

X_Fit  =  Rnd_113 ( : , 5) ; 

V_Fit  =  Rnd_113 ( : , 6) ; 

%. . . .  Brake  Geometry  Data 
Axi  =  BL_46569 ( : , 1) ; 

Dia  =  BL_46569 ( :  ,  2)  ; 

AreaB  =  BL_46569 ( :  ,  3)  ; 

o, 

o 

%....  Set  constant  values  for  brake  geometry  and  fluid  properties 

O, 

o 

gg  =  386.4; 

P_test  =  pr; 
wH20  =  62.4/1728; 
w  =  0 . 853*wH20; 
nu  =  .000111646; 
mu  =  nu* (gg/w) ; 

Cyl_ID  =  3.312; 

ThS_0D  =  1.748; 

AreaA  =  (pi/4 . ) * (Cyl_IDA2  -  ThS_0DA2); 

Orif  =  1.13; 

Disc  =  0.70; 

o, 

o 

%....  Calculate  variable  orifice  area 

o, 

o 

AreaO  =  (pi/4 .)* (Orif A2  -  Dia.*Dia); 
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Wperm  =  pi* (Orif  +  Dia) ; 
Hyd  Rad  =  4*AreaO . /Wperm; 


O, 

o 

%....  Calculate  poly  fits  for  data 

Q, 

O 

Cpl  =  polyf it (t, P_test, Np) ; 

P_Fit  =  polyval (Cpl, t) ; 

o, 

o 

%....  Interpolate  all  calculations  to  travel 

o, 

o 

PX_F i t  =  P_test; 

VX_F i t  =  V_Fit; 

AO_Fit  =  interpl (Axi , AreaO, X  Fit); 

AB  Fit  =  interpl (Axi , AreaB,  X  Fit); 

Hyd  R  =  interpl (Axi , Hyd  Rad,X  Fit); 

O, 

o 

%....  Determine  flow  rate  from  Annular  Chamber 

O, 

o 

Qdot  =  AreaA. *VX_Fit; 

O, 

o 

%....  Calculate  terms  used  in  equation  #4 

o, 

o 

T1  =  2*gg/w; 

T2  =  (Disc. *AreaA. *AB_Fit) .* (Disc. *AreaA.*AB_Fit) ; 

T3  =  AreaA. *AreaA  -  (Disc . *AB_Fit) .* (Disc . *AB_Fit) ; 

QdotB  =  sqrt ( (T2 . /T3) . *T1 . *PX_Fit) ; 

O, 

o 

%....  Calculate  terms  used  in  equation  #6 

o, 

o 

QdotO  =  Qdot  -  QdotB; 

O, 

o 

%....  Calculate  terms  used  in  equation  #5 

o, 

o 

AreaR  =  AreaA. *( 1 . /AO_Fit) ; 

V2  =  AreaA*VX_Fit; 

Reyn  =  (1/mu) * (V2 . *Hyd  R)  ; 

T5  =  AreaAA2*Tl *PX_Fit . / (QdotO . * QdotO) ; 

CdO  =  real (AreaR. / (sqrt (T5  +  1))); 

o, 

o 

%....  Reject  CdO  values  above  1.0 

o, 

o 

%for  j  =  1: size (CdO) 

%  if  CdO(j)  >1.0 

%  CdO(j)  =  1.0; 

%  end 

%end 

ho  =  figure; 

plot (t, x, ' b : ' , t, X_Fit, ' g- ' ) 

%axis (  [0  120  0  25] ) 

title  ({  'XM360  RECOIL  TRAVEL  DATA  &  FIT','ATD-4  TEST  ROUND:  104  '  ,  ' AMSRD-AAR- 

AEW-C ' ,  ' JANUARY  2007  '  } ) 

xlabel('TIME  (m-secs) ') 

ylabel (' TRAVEL  (inches) ') 

grid 

plottools 

hold 
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JANUARY  2007'}) 


hi  =  figure; 


plot (t, P_test, ' b- ' , t, P_Fit, ' g- ' ) 
title ({ 'XM360  BRAKE  PRESSURE  DATA 
xlabel('TIME  (m-secs) ') 
ylabel ( ' PRESSURE  (psi) ' ) 
grid 

plottools 

hold 

o, 

o 

o, 

o 

o, 

o 

h2  =  figure; 
plot (X_Fit, CdO) 

%axis ( [5  20  .8  1.000] ) 

title ({ 'XM360  CALCULATED  VALUE  FOR 

xlabel (' TRAVEL  (inches) ') 

ylabel (' DISCHARGE  COEFFICIENT  ( — ) 

grid 

plottools 

hold 

O, 

o 

o, 

o 

o, 

o 

h3  =  figure; 
plot (V_Fit, CdO) 

%axis ( [100  450  .8  1.000] ) 

title ({ 'XM360  CALCULATED  VALUE  FOR 

C' ,  ' JANUARY  2007  '  } ) 

xlabel ( 'VELOCITY  (in/sec)') 

ylabel (' DISCHARGE  COEFFICIENT  ( — ) 

grid 

plottools 

hold 

O, 

o 

o, 

o 

o, 

o 

h4  =  figure; 
plot (AreaR, CdO) 

%axis ( [0  120  0  6000] ) 

title ({ 'XM360  CALCULATED  VALUE  FOR 

xlabel ( 'AREA  RATIO  ( — )') 

ylabel (' DISCHARGE  COEFFICIENT  ( — ) 

grid 

plottools 

hold 

O, 

o 

o, 

o 

o, 

o 

h5  =  figure; 
plot (t, V_Fit) 

%axis ( [ 0  1  0  600]) 

title ({ 'XM360  RECOIL  VELOCITY  FROM 
xlabel ('TIME  (m-secs) ') 
ylabel (' VELOCITY  (in/sec)') 


FIT  '  ,  ' 


Cd  from  BRAKE  DATA ',' JANUARY  2007'}) 

) 

Cd  from  BRAKE  DATA ' , ' AMSRD-AAR-AEW- 

) 

Cd' ,  ' JANUARY  2007  '  } ) 

) 


DERIVED  FIT ',' JANUARY  2007'}) 
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grid 

plottools 

hold 

o, 

o 

o, 

o 

o, 

o 

h6  =  figure; 
plot (Reyn, CdO) 

%axis ( [0  120  0  6000] ) 

title ({ 'XM360  CALCULATED  VALUE  FOR  Cd JANUARY  2007'}) 
xlabel ( 'REYNOLDS  NUMBER  ( — )') 
ylabel (' DISCHARGE  COEFFICIENT  ( — )') 
grid 

plottools 


O/P  data  for  SP  plotting 

=  [t, x, X_Fit, P_test, P_Fit, V_Fit, AreaR, CdO] ; 
dlmwrite ( f ilenamel , A  Out, ' delimiter ' , ' \t ' , ' precision ' , ' % . 6f ' ) 


hold 

O, 

o 

o, 

o 

o, 

o 


o, 

o 

o, 

o  .  .  .  . 

o, 

o 

A  Out 
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FORTRAN®  Subroutine  of  the  Variable  Cd  Model 

C....Cdfind:  This  routine  is  used  to  determine  the  dynamic  discharge 
C  coefficient  thru  the  maim  orifice  of  the  Schneider-Type 

C  recoil  brake.  It  is  emperical  in  nature  in  that  it  uses 

C  recults  from  a  series  of  steady  state  CFDesign  inns  in 

C  which  the  flow  speed  and  area  ratio  was  varied. 

C 

C...  ORIGINAL  DATE:  01/01/2007 
C...  PROGRAMMER:  R.  G.  GAST 
C 

C....  COMMON  BLOCKS:  None 

C 

C 

C....  INPUT  VARIABLES: 

C  Qdot:  Volumetric  flow  rate  (inA3/s) 

C  Arat:  Area  ratio  of  flow  resevoir  to  restrictive  orifice  (--) 

C 

C  Cd:  Velocity  discharge  coeficient  used  in  Bernoulli  eqn  (--)) 

C 

C 

C....  OTHER  ROUTINES  CALLED:  NONE 
C 

C . 

c 

SUBROUTINE  Cdfind(Qdin,Arin,Cd) 

C 

C....  SET  SPECIFICTIONS,  DEFINE  COMMON  AREAS,  DIMENSION  MATRICES 
C 

implicit  real*  8  (A-H.O-Z) 

dimension  QQ(100),AA(  1 00),CC(  100, 1 00) 

DATA  ktrip  III 
c 

c....  initial  excursion  read  data 
c 

If(ktrip  .eq.  1)  then 
c 

c....  Read  area  ratio  values  and  save  as  common  logarithm 
c 

do  15  Ka=  1,100 
read(28,*)Aval 
if  (Aval  .GE.  0.0)  then 
AA(Ka)  =  log  10(  Aval) 
else 

KA_fin  =  Ka  -  1 
go  to  1 7 
end  if 

15  continue 

1 7  continue 

c 

c....  Read  flow  rate  values 
c 

do  10  Kq=  1,100 
read(28,*)Qval 
if  (Qval  .GE.  0.0)  then 
QQ(Kq)  =  Qval 
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else 

KQ  Ini  =  Kq  -  1 
go  to  12 
end  if 

10  continue 

12  continue 

c 

c....  Read  discahrge  coefficient  values 
c 


do  20  Ka  =  l,KA_fin 
do  20  Kq  =  1  ,KQ_fin 
read(28,*)CC(Ka,Kq) 

20  continue 
ktrip  =  2 
end  if 
c 

c....  Reset  to  local  variables  and  adjust  to  be  within  bounds  of  data 
c 

Arat  =  loglO(Arin) 

Qdot  =  Qdin 
c 

if(Arat  .LT.  AA(1))  then 
Arat  =  AA(1) 

else  if  (Arat  .GE.  AA(KA_fin))  then 
Arat  =  AA(KA_fin)  -,0001*(AA(KA_fin)  -  AA(KA_fin-l)) 
end  if 

if(Qdot  .LT.  QQ(  1 ))  then 
Qdot  =  QQ(l) 

else  if  (Qdot  .GE.  QQ(KQfm))  then 
Qdot  =  QQ(KQ  fin)  -,0001*(QQ(KQ_fin)  -  QQ(KQ  fin-l)) 
end  if 
C 

c....  Search  algorithm  for  Q  &  A  bounds 
c 

do  30  ka  =  1,  KA  fin  -  1 

if(Arat  .GE.  AA(ka)  .AND.  Arat  .LT.  AA(ka+l))then 
KAlow  =  ka 
A1  =  AA(ka) 

A2  =  AA(ka+l) 
go  to  35 
end  if 

30  continue 

35  continue 

do  40  kq  =  1,  KQ  fin  -  1 

if(Qdot  .GE.  QQ(kq)  .AND.  Qdot  .LT.  QQ(kq+l))then 
KQlow  =  kq 
Q1  =QQ(kq) 

Q2  =  QQ(kq+l) 
go  to  45 
end  if 

40  continue 

45  continue 

c 

c....  Linear  search  algorithm  for  Cd  value 
c 
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delA  =  A2  -  A1 

W1  =  sqrt(((Arat-Al)/delA)**2  +  ((Qdot-Ql)/delQ)**2) 

W2  =  sqrt(((Arat-A2)/delA)**2  +  ((Qdot-Q2)/delQ)**2) 

Wtot  =  W1  +  W2 

W22  =  Wl/Wtot 

Wll  =W2/Wtot 

Cl  1  =  CC(KA_low,KQ_low) 

C12  =  CC  (KAlo  w,KQ_lo  w+ 1 ) 

C21  =  CC(KA_low+l,KQ_low) 

C22  =  CC  (KA  lo  w+ 1  ,KQ_low+ 1 ) 


dC  1  dQ  =  (C 1 2-C 1 1  )/delQ 
dC2dQ  =  (C22-C21)/delQ 

dCldA  =  (C21-C1  l)/delA 
dC2dA  =  (C22-C12)/delA 

Cdl  =  Cl  1  +  dC  1  dA*(Arat- A 1 )  +  dCldQ*(Qdot-Ql) 
Cd2  =  C22  +  dC2dA*(Arat-A2)  +  dC2dQ*(Qdot-Q2) 
Cd  =  Wll*Cdl  +  W22*Cd2 


RETURN 


